What is wrong with merely 28% tophat2 mapped reads are counted by HTSeq-count

[bbl]

Here is how my data has been processed: 150bp PE RNAseq, trimmed, tophat2 mapped on both transcriptome and genome using default setting apart from insert-size and std (calculated from sampling), then HTSeq-count using:

htseq-count -q -a 5 --mode=union --stranded=no --type=exon --
idattr=gene_id\
$indir/$sample/accepted_nsort.sam \
$ref/Homo_sapiens.GRCh37.72.gtf > $outdir/$sample.htseq.counts

What worries me is there is only around 28% of the total unique mapped from tophat2 seem to have been properly counted by htseq-count. Could any one divulge what going on here?

Summary of unique mapped reads from tophat2 using samtools flagstat is:
60657139 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
60657139 + 0 mapped (100.00%:-nan%)
60657139 + 0 paired in sequencing
31424549 + 0 read1
29232590 + 0 read2
53208062 + 0 properly paired (87.72%:-nan%)
55095536 + 0 with itself and mate mapped
5561603 + 0 singletons (9.17%:-nan%)
365348 + 0 with mate mapped to a different chr
71148 + 0 with mate mapped to a different chr (mapQ>=5)

Bottom of htseq-count output:
no_feature 11321706
ambiguous 1434804
too_low_aQual 0
not_aligned 0
alignment_not_unique 13865867

The sum of all reads counts by gen_id by htseq-count using awk is 17,068,532.

[kmcarr]

…..
What worries me is there is only around 28% of the total unique mapped from tophat2 seem to have been properly counted by htseq-count. Could any one divulge what going on here?

Summary of unique mapped reads from tophat2 using samtools flagstat is:
60657139 + 0 in total (QC-passed reads + QC-failed reads)
….
Bottom of htseq-count output:
no_feature 11321706
ambiguous 1434804
too_low_aQual 0
not_aligned 0
alignment_not_unique 13865867

The sum of all reads counts by gen_id by htseq-count using awk is 17,068,532.

You are being confused by the way TopHat reports alignments and htseq reports counts. TopHat reports the number of reads aligned, meaning it counts each read of a pair individually so ~60M reads. htseq counts the number of fragments aligned to each gene, meaning it counts each read pair just once. So htseq is basing its counts on ~30M fragments; 17M fragments aligned uniquely would work out to ~57% unique alignment rate.

[bbl]

Thanks kmcarr- my ignorance, didnt remember htseq counts the pair-read only once. Nevertheless, is it reasonable to have 57% of mapped read-pairs for DE analysis? It seems to me still quite a lot loss.

[Brian Bushnell]

That depends on your genome completeness, read quality, and presence of contaminants. Did you do any quality-control? Some basic steps like adapter-trimming, quality-trimming, removal of common contaminants/spike-ins (phiX, etc) can greatly improve the mapping rate. Also, BBMap will map a higher percentage of reads than Tophat, especially if the data is low quality.