subread aligner problem

[pedrocrisp]

Hi, I aligned my reads using subread and I get the following error when running ngs.plot.r. The output files are created but they contain no data, any thoughts?

Loading R libraries.....Done
Analyze bam files and calculate coverageWarning message:  
In headerIndexBam(bam.list) :
  Aligner for: Sample_277_1_WT/Sample_277_1_WT.bam cannot be determined. Style of
standard SAM mapping score will be used. Would you mind submitting an issue
report to us on Github? This will benefit people using the same aligner.
Warning message:
'isNotPrimaryRead' is deprecated.
Use 'isSecondaryAlignment' instead.
See help("Deprecated")
........................................................................................................................................
.............................................Done
Plotting figures...Done
Saving results...Done
Wrapping results up...Done
All done. Cheers!

Thanks peter

[jdblischak]

I was able to get output using BAM files created with Subjunc, which is another short read mapper distributed with Subread.

@pedrocrisp, I could potentially help. Could you please provide some more information?

• What was the command that you ran?
• What is the distribution of quality scores in your BAM file? ngsplot by default removes reads with quality less than 20. The lowest score in my file is 48, so I doubt this is the problem. Here's how to quickly check: samtools view filename.bam | cut -f5 | sort | uniq -c.
• What are the names of your chromosomes? If there is discordance between the chromosome names in your BAM file and the ngsplot annotation data, then there will be no enrichment. You can check quickly with samtools view -H filename.bam.

---

This second part is more for @lishen.

I also receive the message about the mapping scores. Subread has a unique way of calculating the mapping scores. From the User's Guide: The mapping quality score "is a read-length normalized value and it is in the range [0, 60)."

I see that ngsplot using the @PG header to determine the mapper used and thus how to interpret the quality scores (code). The @PG header when mapping with Subread or Subjunc looks like this:

@PG ID:subread  PN:subread  VN:1.5.0-p1

Please let me know if there is anything I can do to help implement support for Subread.

[pedrocrisp]

Hi John, Thanks for the reply. Yes I can see this being a chromosome names issue. Incidentally, I also used subjunc rather than subread.

I ran:

ngs.plot.r -G Tair10 -R genebody -C Sample_277_1_WT.bam -O RRGS_test_genebody -D ensembl -F rnaseq

Quality is good, lots of reads over 20:

samtools view filename.bam | cut -f5 | sort | uniq -c

>
47022898 59
    925 6
    693 7
      1 78
    287 8
   1124 9

Chromosome names do not match:

samtools view -H filename.bam

>
@HD     VN:1.0  SO:coordinate
@SQ     SN:Chr1 LN:30427671
@SQ     SN:Chr2 LN:19698289
@SQ     SN:Chr3 LN:23459830
@SQ     SN:Chr4 LN:18585056
@SQ     SN:Chr5 LN:26975502
@SQ     SN:ChrC LN:154478
@SQ     SN:ChrM LN:366924
@PG     ID:subread      PN:subread      VN:1.4.6

ngsplotdb.py chrnames Tair10 ensembl
>
chr1
chr2
chr3
chr4
chr5
chrM
chrPt

Do you know if there is a work around? Can an alias file be set up for chromosome names or do I need to rename the chromomsomes in all my bam files?

Cheers, peter

[jdblischak]

I'd recommend first renaming the chromosomes in just one of your BAM files and then running ngsplot. If this works, then you'll know that you've isolated the problem. You should be able to do this quickly using samtools (ex):

samtools view -H filename.bam | sed -e 's/Chr/chr/g' | sed -e 's/chrC/chrPt/' | samtools reheader - filename.bam > filename-new.bam
samtools index filename-new.bam

If this is indeed the only problem, and you don't want to rename the chromosomes of all your BAM files (I don't blame you; I also didn't do this), then the only option is to edit the source code.

• Comment out the code inside the function chrTag and have it always return FALSE.
• Change line 659 in the function doCov to replace "chr" with "Chr": v.chrom <- sub('chr', 'Chr', v.chrom)

For background on what this is changing, see Issue #59 that I created and the description of the edits I made to the source code to work with my BAM files. I wasn't able to explicitly test this recommendation since I don't have access to an example file. Give this a try, but if it doesn't work you'll need to share an example BAM for testing for me to help further. Good luck!

[pedrocrisp]

Thanks John, Renaming the chromosomes in my bam, worked a treat, fixed the issue.

I havent tried editing the source yet, I'm testing ngsplot out a little more to see if it suits our needs first.

Cheers