How does HTSeq-count deal with transposons/pseudogenes

[flobpf]

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:

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

[chadn737]

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

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.

[flobpf]

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

[Simon Anders]

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:

I had used the following command-line:

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.

[flobpf]

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!

[chadn737]

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.