(Ported from stack overflow)

I'm trying to summarize gene count using htseq-count; and it's returning 0 counts at every gene. I'm not sure what I'm doing wrong; I think it has to do with the gene flag I'm using. I've tried using the GTF for Arabidopsis TAIR 10: which I got from: ftp://ftp.ensemblgenomes.org/pub/release-46/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gtf.gz as well as the GFF file available on Araport.

I'm using various BAM files I know to have gene expression values.

I've tried both locally and on galaxy The local commands I've tried so far are:

python3 -m HTSeq.scripts.count --idattr=gene_id #bam #gtf > summary.txt
python3 -m HTSeq.scripts.count -i gene_id #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i locus_tag #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i ID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i SeqID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i Parent #bam #gff > summary.txt

I've done the same on galaxy; inputing the -i value in.

Any help as to where I'm going wrong would be great.


  • $\begingroup$ hey @sdshinghal, check a few things. 1. ensure bam file is ok (you can look at it), 2. like the answer below, check your gtf matches the reference genome $\endgroup$
    – StupidWolf
    Mar 6 '20 at 14:09

Maybe the contig names in your BAM files don't match those in the GFF file? Ensembl uses contig names like 1, some other reference sources use chr1.

Generally, for ensembl GFF files, HTSeq's default value for idattr should be fine though.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.