SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
HTSeq-count large proportion of mapped reads with no_feature edm1 Bioinformatics 4 05-26-2014 09:28 AM
Tophat2 Bowtie2 Htseq-count for bacteria chickenmcfu Bioinformatics 2 10-16-2013 06:31 AM
[HTSeq-Count] How are reads counted? syintel87 Bioinformatics 0 06-17-2013 02:18 PM
Discrepancy between HTSeq-count counts and total mapped reads M_staats Bioinformatics 1 03-21-2013 07:16 AM
htseq count missing last mapped fragment mapardo RNA Sequencing 2 11-08-2011 06:27 AM

Reply
 
Thread Tools
Old 06-19-2014, 03:17 AM   #1
bbl
Member
 
Location: London

Join Date: Jul 2010
Posts: 16
Question What is wrong with merely 28% tophat2 mapped reads are counted by HTSeq-count

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.
bbl is offline   Reply With Quote
Old 06-19-2014, 05:11 AM   #2
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,136
Default

Quote:
Originally Posted by bbl View Post
..
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.
kmcarr is offline   Reply With Quote
Old 06-19-2014, 06:57 AM   #3
bbl
Member
 
Location: London

Join Date: Jul 2010
Posts: 16
Default

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.
bbl is offline   Reply With Quote
Old 06-19-2014, 11:00 AM   #4
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

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.
Brian Bushnell is offline   Reply With Quote
Reply

Tags
htseq-count, rnaseq, tophat2 mapping

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:13 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO