This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Difference Between Htseq-Count And Bedtools Multicov

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.

                        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

rna-seq

you should show how you used them

2 answers

HTseq-count drops reads that can not be assigned unambigously based on program parameters deciding rules http://www-huber.embl.de/users/anders/HTSeq/doc/count.html.

  • If you have isoforms of a gene with overlapping exons, htseq-count will discard reads that fall into these exons. You should find the count of reads in the 'ambigous" output row.

  • bedtools will assign them to every isoform where there is an overlap, you should have more counts than reads for the same annotation where you have htseq-ambigous reads.

HTseq-count is said to be able to assign reads to ensembl-GTF file genes by aggregating over isoforms combined by the "gene_id". I was unlucky until now to make this work well.

For that, you can try htseq-count xxx.sorted.sam xxx.gtf > gene_counts.txt

I think the python script in DEXSeq can also give the same result

you can check your GFF/GTF file and htseq-count command line:

1. check the col-3 (feature field) of GFF/GTF file, that should contain "exon"

2. in your htseq-count command line: --type=FEATURETYPE, default is "exon"

Log in to answer this question.