Samtools/htseq-count numbers don't add up.

[tirohia]

I have aligned a bunch of paired end samples to a reference genome using tophat. For each sample I then get an accepted_hits.bam and an unmapped.bam file.

When I use

samtools view -c -F 4 accepted_hits.bam

I get 38216419 reads mapped, which is fine. It's what I should expect, all the reads in accepted hits are mapped to the genome. I can drop the -F 4 and get the same number.

When I use htseq-count though:

htseq-count -f bam -r name -o codingRNA -t gene -i ID accepted_hits.bam actinidia.gff

I get 21593773 (by recording the output and adding up all the counts). This includes __no_feature, __ambiguous, etc. Why the two different numbers? It's not as if one of them if individual reads and the other read total, as one isn't twice the other. Surely if everything is mapped to the reference genome, then the total from htseq-count should be the same as the samtools count, given that the htseq-count should include all the counts mapped to non-coding regions in the genome?

[dpryan]

Given that you used tophat2, it's likely that this is due to a few different reasons. Firstly, tophat2 will produce singletons, which will be counted only once by both samtools and htseq-count but would throw off any "half the number of reads should be the sum of the counts" rule of thumb. Secondly, it's quite likely that you have multiple alignments from some reads in the BAM file (I don't recall what tophat2 defaults to anymore). Have a look at samtools idxstats and that should clarify things.

[tirohia]

Thanks for the reply. It does give me another question or two though.

The bit about singletons makes sense. Does that raise any problem with raw read counts? If I have some of the reads with only one end mapped and some with both I mean?

Do you mean I might have reads in the bam file mapping to multiple places? I would have thought htseq-count would have listed those as ambiguous or some such.

As far as samtools idxstats goes, what am I looking for in? I understand that it's meant to be a list of Chr name, length, reads mapped and unmapped. Though I don't know how that would work give unmapped reads are genome wide rather than chromosome level. Unless reads mapped to other chromosomes are meant to be counts as unmapped to the chromosome you're looking at.

I don't see how it would let me ascertain if there have been reads that have been aligned to multiple places.

Either way, mine is rather odd.

Chr1 17971566 598636 0
Chr10 15490421 426329 0
Chr11 17082262 1398852 0
Chr12 13002763 4185681 0
Chr13 17372821 624540 0
Chr14 12031045 296495 0
Chr15 19001909 682029 0
Chr16 11476499 437375 0
Chr17 14675105 715411 0
Chr18 14814751 997682 0
Chr19 11145676 2113419 0
Chr2 15003700 515401 0
Chr20 14375345 407763 0
Chr21 14287830 476878 0
Chr22 12982530 383782 0
Chr23 20688064 926548 0
Chr24 17702303 595310 0
Chr25 14321294 1934664 0
Chr26 17153587 621163 0
Chr27 10917718 299124 0
Chr28 16406049 836221 0
Chr29 14180785 478150 0
Chr3 20519224 1308238 0
Chr4 13065296 2270267 0
Chr5 18825850 1363623 0
Chr6 18208465 883421 0
Chr7 18704015 1445303 0
Chr8 21083321 693624 0
Chr9 11598160 921655 0
Unknow 177361715 9378835 0

I'm presuming that this means the accepted_hits.bam file has nothing but mapped genes in it?
Which would make sense, given that the unmapped reads are in a different file. It doesn't clarify ,much for me though. Apart from the fact that it all adds up to the same number that samtools view -c gives me.

[dpryan]

Singletons won't adversely affect htseq-count.

I actually meant flagstat, though I typed idxstats, mea culpa. You want to look at the number that aligned as singletons and similar metrics.