NH:n:x in BAM output

[CPCantalapiedra]

Maybe my last sentence was not right after all. It seems that NH:i:1 and MAPQ50 are constant when changing -g (except if -g, in which case all the reads are (erroneously) reported as NH:i:1 and MAPQ50). However, NH:i (>1) and MAPQ (<50) counts do not coincide with align_summary multiple mapping reads number. I would really like to understand align_summary report.

[CPCantalapiedra]

Dear all,

... Tophat2 reports the NH flawlessly, therefore there must be a way to run Bowtie(2) with an option to generate this report.

Federico Giorgi (CUMC) & Davide Scaglione (IGA)

I am not sure it reports it flawlessly. When tuning the -g (max-multihits) option, the number of NH:i:>1 changes, and this seems wrong to me.

From the TopHat2 manual:
"If there are more alignments with the same score than this number, TopHat will randomly report only this many alignments"

But not reporting them should not mean that they are not there, as align_summary.txt seems to understand. However, in my case (TopHat2.0.9) all reads with multiple mappings > -g are reported with MAPQ50 and NH:i:1

[wks2105]

I agree that this problem needs addressing ..

[giorgifm]

To remove multiple hit reads, I am also using
samtools view -q4 -F4 file.bam

Following the explanation given by this post:

http://seqanswers.com/forums/archive/index.php/t-19491.html

I was aware of the solution by Simon Anders, and I pray the gods every day for its existence. However, -F4 -q4 may remove primary, unambiguous but low quality mappings, and therefore this is an efficient but patchy solution

[giorgifm]

while XS:i the score for the second-best one.

Dear Gianza, so if an alignment row contains the XS field, the read is by definition aligning in more than one place?

Or rather, if the reported AS score is equal to the Q field of the BAM, then we are looking at (one of) the best hits? In this case the problem remains, and it would be simpler to just count the number of reads/fragments and check for their Q scores.

So what I believe now is that the only way to obtain the NH is writing an addNH script, that processes the entire readname-sorted BAM to add the NH field based ont he Qs and the FLAGs contained within the BAM itseld

[peer.b]

To remove multiple hit reads, I am also using
samtools view -q4 -F4 file.bam

Following the explanation given by this post:

http://seqanswers.com/forums/archive/index.php/t-19491.html

I concur that the lack of an NH:i:n field in BWA, and especially in Bowtie is unbelievable, I see no reason why the programmers omitted this field. If it weren't for "standards" and complaining reviewers, I would switch to other aligners like ERNE at once

[Gianza]

Reading the Bowtie2 manual, it seems that the AS:i: flag reports the score for the best alignment, while XS:i the score for the second-best one. Thus one could compare these two value to figure out wheather a read is mapping ambiguosly or not.

Still, I'm really frustated not being able to control the maximum number of mismatches and/or gaps as tophat2 easily does with two dedicated options.

Cheers

Davide

[giorgifm]

Dear all,

the option specified for the number of hits in a SAM output, NH, is amazingly missing from the two most popular aligners (namely BWA and Bowtie). However, Tophat2 reports the NH flawlessly, therefore there must be a way to run Bowtie(2) with an option to generate this report.

So far, too much of the community is resorting to the samtools view -q3 trick, which is able to extract alignments with a mapping quality below 50% (by definition, a multiple hitting read will have all its n equally good alignments with a score of at most 1/n).

Thanks for any hints. I hope this post will become viral.

Federico Giorgi (CUMC) & Davide Scaglione (IGA)