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