SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
htseq-count for sam and gff3 sofia17 RNA Sequencing 45 11-04-2016 03:32 PM
Strange error when using htseq-count shhuang Bioinformatics 13 11-19-2012 12:40 AM
htseq count missing last mapped fragment mapardo RNA Sequencing 2 11-08-2011 05:27 AM
How does HTSeq-count deal with transposons/pseudogenes flobpf Bioinformatics 5 08-17-2011 12:53 PM
Count difference by htseq-cout and samtools DZhang Bioinformatics 2 07-03-2011 11:05 AM

Reply
 
Thread Tools
Old 03-01-2012, 06:48 AM   #1
paolo.kunder
Member
 
Location: Milano, Italy

Join Date: Aug 2011
Posts: 92
Default htseq-count

Hi,
I was wondering about the output of htseq-count to count reads associated with transcripts.

I run htseq-count to estimate reads associated with genes with the command below:

htseq-count --stranded=yes --mode=intersection-nonempty file.sam file.gtf >> gene_id.txt

example I found one gene ENSGxxx with associated 700 reads count.

To check how many reads are associated with the different transcripts of ENSGxxx gene i run the command below:

htseq-count --stranded=yes --idattr=transcript_id --mode=intersection-nonempty file.sam file.gtf >> transcript id.txt

I would expect that the sum of the reads associated with each transcript is almost similar to the reads associated with the gene, but it seems not the case..

example:
gene 700 counts
transcript1 0 counts
transcript2 4 counts
transcript3 72 counts


any suggestions?
Ciao!
paolo
paolo.kunder is offline   Reply With Quote
Old 03-01-2012, 07:05 AM   #2
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

I can't say for certain, but this is what I bet is going on.

Even with the interesection_nonempty option, reads that map equally well to more than will attribute will not be counted and instead will be counted in the ambiguous category.

The different transcripts are likely going to share the majority of the same exons. So by using transcript_id as the attribute, most of those reads are going to be counted as ambiguous since they could be counted equally for any one of the three transcripts. The counts that you do get for each transcript are those that map to the exons unique to each of those transcripts.

Look at the end of your count list at the category "ambiguous" I almost guarantee you that the number of counts in this category is vastly more when you use transcript_id versus gene_id.
chadn737 is offline   Reply With Quote
Old 03-01-2012, 07:12 AM   #3
paolo.kunder
Member
 
Location: Milano, Italy

Join Date: Aug 2011
Posts: 92
Default

exactly
Ambiguos reads in transcript id is much more bigger than gene id,

I tried also to use --mode=union but nothing changed,

so I think is useless to counts reads associated with transcript with htseq-count,

any other program could I use?

Paolo

Last edited by paolo.kunder; 03-01-2012 at 07:16 AM.
paolo.kunder is offline   Reply With Quote
Old 03-01-2012, 07:45 AM   #4
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

What do you want to do with it after you have your counts? Is there a specific program you want to use?

There is RSEM

http://deweylab.biostat.wisc.edu/rsem/
chadn737 is offline   Reply With Quote
Old 03-01-2012, 08:40 AM   #5
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

htseq-count is designed to count for genes, not for transcripts. Post #2 correctly explained why the counts are so low.

Why would you want to count reads per transcript? Given that the vast majority of reads will overlap with more than one transcript, such count values have little to do with transcript expression strengths. (If it were as easy, nobody would bother with developing algorithms for isoform deconvolution.)
Simon Anders is offline   Reply With Quote
Old 03-01-2012, 11:39 PM   #6
paolo.kunder
Member
 
Location: Milano, Italy

Join Date: Aug 2011
Posts: 92
Default

Thanks for the replies, I'll have a look to RSEM and continue my analysis with DEXSeq and so on...




Paolo
paolo.kunder is offline   Reply With Quote
Old 02-26-2014, 09:47 AM   #7
geneart
Member
 
Location: DC area

Join Date: Sep 2011
Posts: 42
Default

Hello Folks,
I had a question about HTSeq output which is kind of in tune with this blog....I have sam files as output as well as the count data in text format. So I wanted to look at the no_feature tags that HTSeq puts out to then match with GFF file to see what those no_feature are ( if they were transposons or intronic regions etc).
My problem is that according to HTSeq count text it generates , i have a value of 33190 for no_feature. However when I either grep for 'no_feature' or reverse grep for 'ENSBTAG' from my sam file to get counts, I get 33055. Why this variation? I would love to get the sequences associated with this no_feature out of sam file to further process it . Any suggestions?
Thanks
Geneart.
geneart is offline   Reply With Quote
Old 02-26-2014, 10:45 AM   #8
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Coincidentally, I spotted a minor bug today when preparing the new release of HTSeq: If a read maps to a chromosome or contig for which no exons at all are listed in the GTF file, the read will be correctly counted as "no_feature" but will not be written to the "samout" file. This is fixed in version 0.6.0, available since an hour or so. Maybe this solves your issue.
Simon Anders is offline   Reply With Quote
Old 02-26-2014, 11:04 AM   #9
geneart
Member
 
Location: DC area

Join Date: Sep 2011
Posts: 42
Default

Thanks very much!! Will rerun the process with new version.
geneart is offline   Reply With Quote
Old 03-11-2014, 12:02 PM   #10
wwq413
Junior Member
 
Location: usa

Join Date: Dec 2010
Posts: 6
Default

Quote:
Originally Posted by Simon Anders View Post
Coincidentally, I spotted a minor bug today when preparing the new release of HTSeq: If a read maps to a chromosome or contig for which no exons at all are listed in the GTF file, the read will be correctly counted as "no_feature" but will not be written to the "samout" file. This is fixed in version 0.6.0, available since an hour or so. Maybe this solves your issue.
You are so right that it worked after i changed into exon in .gff. Thanks!
__________________
wwq
wwq413 is offline   Reply With Quote
Old 10-22-2014, 04:45 AM   #11
Gonza
Member
 
Location: Ithaca, NY

Join Date: Mar 2013
Posts: 78
Default

Dear All,

I am trying to run htseq-count but it keep failing.

Script:
htseq-count --stranded=no -o samout_test accepted_hits.sam Arabidopsis_thaliana.TAIR10.23.gtf

100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
485101 GFF lines processed.
Error occured when processing SAM input (line 20 of file accepted_hits.sam):
('need more than 10 values to unpack', 'line 20 of file accepted_hits.sam')
[Exception type: ValueError, raised in _HTSeq.pyx:1277]

I converted the tophat2 file accepted_hits.bam to accepted_hits.sam and used this as input. I am really struggling to understand what this means and how to solve the problem. Any ideas?

thanks in advance
G
Gonza is offline   Reply With Quote
Reply

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 01:12 PM.


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