htseq-count

[paolo.kunder]

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

[chadn737]

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.

[paolo.kunder]

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

[chadn737]

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/

[Simon Anders]

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.)

[paolo.kunder]

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




Paolo

[geneart]

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.

[Simon Anders]

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.

[geneart]

Thanks very much!! Will rerun the process with new version.

[wwq413]

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!

[Gonza]

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