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 04:32 PM
Error with GTF file when using htseq-count MDonlin Bioinformatics 13 01-13-2015 09:29 AM
Strange error when using htseq-count shhuang Bioinformatics 13 11-19-2012 01:40 AM
htseq count missing last mapped fragment mapardo RNA Sequencing 2 11-08-2011 06:27 AM
Count difference by htseq-cout and samtools DZhang Bioinformatics 2 07-03-2011 12:05 PM

Reply
 
Thread Tools
Old 08-15-2011, 01:32 PM   #1
flobpf
Member
 
Location: USA

Join Date: Apr 2010
Posts: 76
Default How does HTSeq-count deal with transposons/pseudogenes

Hi,

I noticed that the HTSeq-count output does not have any reads unambiguously mapped to pseudogenes and transposons. I havent checked the no-feature and ambiguous reads. I wonder whether HTSeq classifies reads mapping to pseudogenes/transposons as no-feature or as ambiguous? If neither, then what?

I had used the following command-line:

Code:
python -m HTSeq.scripts.count -m union -s no -t gene -i ID -o myfile.gff.readcounts accepted_hits.sam myfile.gff
Any suggestions appreciated!

Thanks
Flobpf
flobpf is offline   Reply With Quote
Old 08-15-2011, 02:45 PM   #2
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

I'm sure Simon Anders will probably correct me on this, but I would assume its your GFF file.

My experience is that HT-seq-count is picky about what features it uses and will only use something if it is named as an exon. However, some GFF files use different names for psuedogenes or transposons exons. To give an example, here is a bit from the TAIR10 GFF. Exons for psuedogenes are listed as "pseudogenic_exon."

Quote:
Chr1 TAIR10 pseudogene 402693 402961 . - . ID=AT1G02136;Note=pseudogene;Name=AT1G02136
Chr1 TAIR10 pseudogenic_transcript 402693 402961 . - . ID=AT1G02136.1;Parent=AT1G02136;Name=AT1G02136.1;Index=1
Chr1 TAIR10 pseudogenic_exon 402886 402961 . - . Parent=AT1G02136.1
Chr1 TAIR10 pseudogenic_exon 402693 402808 . - . Parent=AT1G02136.1
Chr1 TAIR10 gene 403190 404456 . - . ID=AT1G02140;Note=protein_coding_gene;Name=AT1G02140
Chr1 TAIR10 mRNA 403190 404456 . - . ID=AT1G02140.1;Parent=AT1G02140;Name=AT1G02140.1;Index=1
Chr1 TAIR10 protein 403467 404401 . - . ID=AT1G02140.1-Protein;Name=AT1G02140.1;Derives_from=AT1G02140.1
Chr1 TAIR10 five_prime_UTR 404402 404456 . - . Parent=AT1G02140.1
Chr1 TAIR10 CDS 404186 404401 . - 0 Parent=AT1G02140.1,AT1G02140.1-Protein;
When I wanted to use HT-seq to count how many reads mapped to one of these, I first had to modify my GFF file, so that "psuedogenic_exon" was listed as "exon."

I would check your GFF file to see if they are annotated in the same way.
chadn737 is offline   Reply With Quote
Old 08-16-2011, 10:10 AM   #3
flobpf
Member
 
Location: USA

Join Date: Apr 2010
Posts: 76
Default Indeed

Hi chadn737,

Thanks for the reply!

I understand that these two features are not defined as "genes". I'm using TAIR10 too and as you rightly mention, these are defined as "Pseudogenes" and "transposable_element_gene". However, what I'd like to know is whether HTSeq classifies reads mapping to these features as "no_feature" or "ambiguous".

UPDATE: I just finished analyzing the no_feature reads. Several no_feature reads were mapping to pseudogenes and transposons. I think it would be better if the program classified these in some other category than no_feature to prevent misunderstanding.

Thanks
flobpf is offline   Reply With Quote
Old 08-17-2011, 02:04 AM   #4
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Quote:
I think it would be better if the program classified these in some other category than no_feature to prevent misunderstanding.
The GFF format has a very loose (and partly even contradictory) definition. Especially, there is no controlled vocabulary for the "type" field in the third column, and hence,
different providers of GFF file may use different terms. Hence, it has to be left to the user to sort this out.


By the way:

Quote:
I had used the following command-line:
Code:
python -m HTSeq.scripts.count -m union -s no -t gene -i ID -o myfile.gff.readcounts accepted_hi
You have explicitly instructed htseq-count to only look at the "gene" features, rather than the "exon" features. I assume you intended to count reads mapping to introns.
Simon Anders is offline   Reply With Quote
Old 08-17-2011, 07:01 AM   #5
flobpf
Member
 
Location: USA

Join Date: Apr 2010
Posts: 76
Default Makes sense...

Thanks Simon. Your logic makes sense.

Also, I did intend to get all reads mapping to introns (actually, I wanted all "intergenic" reads, even excluding those mapping to pseudogenes/transposons...thus the question...).

Thanks!
flobpf is offline   Reply With Quote
Old 08-17-2011, 01:53 PM   #6
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

Quote:
Originally Posted by flobpf View Post
Thanks Simon. Your logic makes sense.

Also, I did intend to get all reads mapping to introns (actually, I wanted all "intergenic" reads, even excluding those mapping to pseudogenes/transposons...thus the question...).

Thanks!
If all you are interested in is the number of reads mapping to intergenic regions, I can help you because I had the same question with some of my data. I ended up creating a GTF file from the TAIR10_intragenic sequence file that could be used by HTSeq-count.
chadn737 is offline   Reply With Quote
Reply

Tags
annotation, gff, htseq-count, illumina

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 08:45 PM.


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