SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
htseq-count low count problem gandalf886 Bioinformatics 3 08-23-2014 08:05 AM
Problem with HTSeq-count anikng RNA Sequencing 3 08-16-2013 09:33 PM
HTSeq-count issue zuco RNA Sequencing 2 08-16-2013 12:30 AM
samtools sorting problem yaximik Bioinformatics 9 04-22-2013 06:22 AM
Issue with htseq-count cpleis Bioinformatics 8 10-15-2012 10:31 AM

Reply
 
Thread Tools
Old 10-28-2014, 04:22 AM   #1
bbl
Member
 
Location: London

Join Date: Jul 2010
Posts: 16
Default samtools sorting issue or HTSeq-count problem?

Hello all,

This has puzzled me for days and I couldnt find any explanation on the internet.
I noticed HTSeq-count gives the reads counts of a gene as 0. But when viewed in IGV with the accepted_hits.bam from tophat2 alignment, there are hundreds to thousands reads aligned in that gene region (from different samples).

This is how I used HTSeq-count from HTSeq/0.6.1p1 and samtools/1.1
samtools sort -T tmp -o accepted_hits_nsort.bam -n accepted_hits.bam
htseq-count --format=bam --strand=no --order=name -a 5 --mode=union accepted_hits_nsort.bam Homo_sapiens.GRCh37.72.gtf > htseq_count_0.6.txt

However, when I break down the accepted_hits.bam and extract the chromosome where the gene is, HTSeq-count gives the right reads counts for this BAM file.

My accepted_hits.bam contains 60,845,216 alginemnt using samtools view -c accepted_hits.bam. I dont think for PE RNA-seq, file at size of 4.5G would be a problem for samtools or HTSeq-count to deal with?

My suspicions are:
1. samtools name sort problem
2. htseq-count bug
3. out of memory

Could any one give any hint here? Any suggestion would be appreciated. Thanks.
bbl is offline   Reply With Quote
Old 10-28-2014, 04:35 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,476
Default

Try using the -o option and see what happens to some of those reads that should be getting counted.

BTW, featureCounts is MUCH faster, which is especially convenient on larger files.
dpryan is offline   Reply With Quote
Old 10-28-2014, 05:32 AM   #3
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,135
Default

What gene are you looking at in IGV? The most likely answer is that you are looking at a duplicated gene and the 100's to 1000's of reads which map to it also map to the other copies elsewhere in the genome. When a read maps to multiple locations HTSeq-count will not count it for any of the genes it maps to, rather it will be be counted among the 'alignment_not_unique'. Your experiment of running HTSeq-count on just one chromosome shows that the additional mappings for these reads fall on other chromosomes. When HTSeq-count is given an incomplete alignment set with only one valid mapping for these reads it will naturally count them for that one alignment. This demonstrates why it is a bad idea to perform read counting on only partial alignment sets, unless you are prepared to deal with these types of situations.

If you want to confirm this identify one of the reads aligned to the gene you have identified in IGV and then search for that read name in your accepted_hits.bam to see if it is multiply mapped.
kmcarr is offline   Reply With Quote
Old 10-28-2014, 08:51 AM   #4
bbl
Member
 
Location: London

Join Date: Jul 2010
Posts: 16
Default

Quote:
Originally Posted by dpryan View Post
Try using the -o option and see what happens to some of those reads that should be getting counted.

BTW, featureCounts is MUCH faster, which is especially convenient on larger files.
Thanks very much dpryan. I shall try featureCounts too.
bbl is offline   Reply With Quote
Old 10-28-2014, 08:57 AM   #5
bbl
Member
 
Location: London

Join Date: Jul 2010
Posts: 16
Default

Thanks very much kmcarr. Indeed what you suspected is what happens!
1st, CEBPA is the gene I'm intereted. From samout by HTSeq-count, reads are mapped to two and three different genes incl. CEBPA. Therefore, they are counted as ambiguous. Actually these reads are also mapped to 1. 'retired' gene, 2. antisense in my .gff ref downloaded from Ensembl.
Since I'm looking at RNAseq samples treated using polyA+tail , I suppose it makes sense to 'manually' modify the .gff ref? For example, remove those pseudogenes, non-coding RNA etc...Would you suggest some way to modify it to be more suitable, if it is a right thing to do?

Many thanks again.
bbl is offline   Reply With Quote
Reply

Tags
htseq-count, reads counts, samtools sort

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 02:47 PM.


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