This is a test version of Biostars. For the public version, visit https://www.biostars.org.
not primary alignment flag filtering does not work

Hello,

I have a bowtie2 result looking like this;

27823422 reads; of these:   27823422 (100.00%) were unpaired; of these:
2464773 (8.86%) aligned 0 times
19083986 (68.59%) aligned exactly 1 time
6274663 (22.55%) aligned >1 times
91.14% overall alignment rate

I am trying to count the uniquely mapped reads, and in this case I think the value should be 19083986 since it should be aligned once.

I have tried these commands;

samtools view -h -F 256 -c input.bam 
samtools view -h -F 0x100 -c input.bam

But they are both giving the total (27823422) number. Why could it be like that?

Thank you in advance.

filtering samtools chipsamples chipseq chipseqanalysis
• 2,130 views • link
updated by jkbonfield jkbonfield • written by buffealo • buffealo

Different tools and people have different notions of what "primary" means. It's not even that clear in the spec. Is a supplementary read part of the primary alignment if it's a second portion of the main read (ie supplementary but not secondary)? I'd say yes, but I think that's a minority view as often people claim primary alignment is non-supplementary too.

I suspect you're mixing up primary vs non-secondary too.

Rather than relying on the summary output from the aligner for diagnosing the issue, you should try samtools flagstat on your file to see what the counts really are. That'll give you an idea which flag values to filter on.

• 0 views • link
written by jkbonfield jkbonfield

1 answer

Count up that total number of reads. It might be that Bowtie is choosing to only report a single alignment for multi-mappers, so there will be no secondary alignments. You might have to use grep or awk to inspect the tags to determine which reads mapped uniquely, instead of the flag.

Even if multi-mappers were being reported multiple times, there would be one instance of each read marked as primary, so filtering for primary alignments only would not remove multi-mapping reads; it would just ensure they are in your output once and only once.

• 0 views • link
written by swbarnes2 swbarnes2

Log in to answer this question.

More posts like this