SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
CuffDiff: genes with highly abundant reads have FPKM 0. uniportdb Bioinformatics 3 02-16-2013 09:09 PM
Gene reads count from RNA-seq sam file by CIGAR Jerry_Zhao Bioinformatics 8 02-01-2013 11:01 AM
count reads or count base-pairs yuelics Introductions 3 07-29-2011 05:41 AM
Quantification: count reads or count base pairs? yuelics Bioinformatics 0 07-27-2011 04:48 AM
filtering out reads from abundant transcripts before using velevt Marco Bioinformatics 2 02-03-2010 09:54 AM

Reply
 
Thread Tools
Old 03-05-2013, 12:20 PM   #1
Babooncanfly
Junior Member
 
Location: USA

Join Date: Mar 2013
Posts: 5
Question hsteq-count giving zero reads on abundant gene

Dear all and Simon,

htseq-count ran fine on my mm9_transcriptome gtf file. No warning messages.

The output .count file looks mostly correct also. However, an abundant gene "Terc" gives zero read.

The gene is expressed, and in ucsc genome browser you can see many reads on the feature (bigwig converted from the same .sam file used in htseq-count). So it's not a mapping issue and anything upstream of .sam is not the source of the problem.

I checked the gtf file and the coordinates of "Terc" matches that of the genome browser so it's not that either.

It's not a strand-thing because I enforced no strand-specificity.

I then isolated the "Terc" line (single exon) and made a mini gtf file, htseq-count against just this one feature, and got ~300 reads, which seems reasonable. But with the rest of the file there it drops to zero.

I then suspected some overlapping feature soaking up the reads that should've gone to Terc, but fives lines above and below the "Terc" feature contain no overlapping coordinates, and the gtf file is sorted by genomic coordinates.

I then checked if the same thing's happening to other genes and at least one other gene I looked at gives the same count if I used the complete gtf file or just itself. This gene does have 4 exons so 4 lines but that's probably not the issue?

At this point I'm at my wits end. I now suspect that reads belonging to Terc in the sam file are also mapped to somewhere else but I a) don't think Terc is homologous completely with some other gene; b) it's hard to imagine that would reduce count from 300 to 0.

Any ideas? I recall seeing someone else having a similar problem of seeing zero counts on select genes that they know are abundantly expressed in their organism (I think it was a bacteria). Any help is appreciated.
Babooncanfly is offline   Reply With Quote
Reply

Tags
error, gtf, htseq, htseq-count

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


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