difference between htseq-count and bedtools multicov

[Simon Anders]

I never really ran into these problems myself and I wonder if this might be because I only use Ensembl annotation. Ensembl changes their gene IDs quite often so that you have to remap these ENBG000... IDs to HGNC symbols anew for each Ensembl release but the advantage is that this gives the Ensembl people the flexibility to change easily what they consider one gene. This might explain why I feel that the described issue with overlapping genes is somehow avoided in Ensembl annotation

[sdriscoll]

maybe it's just me but i've always been annoyed with official gene symbols. from the start i've never thought they were reliable to base any part of my codes around so when i've written analysis tools for my lab i use more checks to make sure my code knows what it's looking at and so i don't confuse anybody. from a programmer's point of view it's just a design constraint and one less assumption that can be made along the way. i end up dealing with a lot of different annotations, sometimes custom made, and different genomes and it's easier for me if my tools don't put too much requirement on the naming conventions of the annotation.

i'm only pointing out the fact that even if you're working with official gene symbols as the 'gene_id' field then a counter like htseq-count will end up throwing out some reads from many loci where all that's actually happening is some illogical historical naming issues. that's not the researcher's fault. some of those locations in the genome are overlap of adjacent genes but sometimes it's just what looks like a typical alternatively spliced gene locus but for some reason not all of the isoforms have the same gene name. those are the ones that are the issue. it's a minor point - yes if you have the right annotation it's irrelevant.

the human genome is even more complex. for example this SMN2 locus shows what I'm talking about.

http://genome.ucsc.edu/cgi-bin/hgTracks?hgHubConnect.destUrl=..%2Fcgi-bin%2FhgTracks&clade=mammal&org=Human&db=hg19&position=chr5%3A70220768-70248842&hgt.positionInput=SMN2+%28Homo+sapiens+survival+of+motor+neuron+2%2C+centromeric+%28SMN2%29%2C+transcript+variant+d%2C+mRNA.%29&hgt.suggestTrack=knownGene&Submit=submit&hgsid=311636999&pix=1330&hgFind.matches=uc003kak.3&knownGene=pack

looks like your typical alternatively spliced locus but it houses SMN1 and SMN2. UCSC's track and the Refseq track both show this naming issue.

[pbluescript]

The UCSC gene name substitution is a problem I've run into before too.
Here's a blog post I found that offers a solution for anyone who's interested:

http://compbionovice.blogspot.com/2012/06/when-trying-to-do-gene-and-transcript.html

[syfo]

Yes, replacing UCSC gene_id by NCBI gene symbols is the good move. Interesting to see that still more than 10% of the historically annotated loci overlap with at least another one in the mouse genome, even when considering "real gene symbols".

[sdriscoll]

That's correct. I however am not working with the raw Gtf to get that number. I swap in the real gene symbols into the gene_id field. I also ran an analysis of the entire gtf collapsing different transcript ids down into loci where different transcripts share some sequence between them. I actually just use gffread (from Cufflinks) for that. Then parsing that result I find >2200 loci with multiple gene names, and obviously multiple transcript ids. If you look at the gene loci I mentioned in the UCSC genome browser for mm9 you can see what i mean.

[syfo]

Sdriscoll, I think that the main issue here is not due to "genuinely" overlapping genes but to the bug Simon mentioned: in a raw GTF from the UCSC table browser all the transcripts have a different "gene_id" whether they come from the same locus or not.
I count 55,419 different identifiers in the "gene_id" field of the UCSC knownGene GTF file for mm9, not 2,200.

[sdriscoll]

Hi Simon - you're right about that. If you use the UCSC GTFs straight away then your program won't work well at all. I was being a little dramatic in my post however even if you have the gene names inserted in your GTF, which I always do, there's still a large number of multi-gene loci where information can be lost because of the unique gene_id requirement in htseq-count. With a quick loci clustering and a search for those with multiple gene names (using UCSC's knowngene annotation for mouse mm9) there's over 2200 loci with more than one gene_id which is close to 10% of the annotated loci in the mouse genome. I don't work with human much but I suspect there's more overlap there. Maybe not all of those multi-gene loci share all of their exons with each other but some of them do. For example there are loci like CC1/Rb1cc1 or Sgk3/Cisk. One would lose quite a few reads in those loci. So a bit more work would have to be done to one's annotation to ensure that they aren't losing too much information with htseq-count. For example htseq-count works beautifully if you collapse these loci into exon-unions. I've written my own read counters over the years and what seems to work best for me is to allow reads to be assigned to multiple features but to weight those reads by 1/(number of features). I could futher restrict that to only allowing those reads to be ones that are shared across features within the same locus but if I'm working with alignments that don't inlcude multi-aligned reads then that doesn't matter. This way I get as many reads assigned to features as possible and when I collapse the loci down into "genes" I can just add the features' counts. In the end I get the same counts as I would have by using a modifed GTF with overlapped features collapsed into unions without having to generate that extra GTF file.

[Simon Anders]

sdriscoll, te source of your problems is simply that you use annotation files that are unsuitable for htseq-count.

If a read alignes two several transcripts of the same gene, it is counted for this gene; only if it overlaps with different genes, it is dropped.

To this end, your GTF file needs to indicate which transcripts belong to the same gene. The GTF standard specifies that exon lines from two transcripts of the same gene should have the same 'gene_id' but different 'transcript_id' values. Unfortunately, the Table Browser function of the UCSC Genome Browser returns faulty GTF files that just copy the transcript ID into the gene ID field, and with these, htseq-count does not work.

So, please don't use the UCSC Table Browser. (They explained on previous occasion that they are unable to fix this bug, which, in my opinion, is rather serious.)

[sdriscoll]

i'd guess it's due to the fact that htseq-count only reports one hit per aligned read assuming that read is aligned uniquely and does not overlap multiple features listed in your GTF. if an aligned read hits more than one feature in your GTF then it doesn't report that hit. bedtools gives you raw hits which includes every 1 hit for every intersection of every alignment with any features in the GTF no matter how many times it aligned or how many features it hit. you might think, "wow, htseq-count is dropping a lot of information". yes, it is! i've moved to using other tools to count hits to genes (RSEM/eXpress) since they disambiguate those ambiguous alignments and as a result you get counts for all of your aligned reads. in a genome with alternative splicing you lose too much data using htseq-count, in my opinion.

[camelbbs]

Hi,

Anyone have used htseq-count and bedtools multicov, I use both them to get the read counts for ref genes. But I got very different read counts between them. Which one is more common..

bedtools htseq-count
uc001aaa.3 8 0
uc001aac.4 2309 0
uc001aae.4 753 0
uc001aah.4 2309 0
uc001aai.1 586 0
uc001aak.3 0 0
uc001aal.1 0 0
uc001aam.4 1315 6