Difference Between Htseq-Count And Bedtools Multicov

Question

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

                        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

[Istvan Albert]

you should show how you used them

Answer 1

[Ido Tamir]

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.

[camelbbs]

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

Answer 2

[wm]

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"