Mapping reads and quantifying genes - Metagenomic workshop

Question

[arshad1292]

Hello,

I am using the following metagenomic workshop tutorial to analyse my own metagenomic data.

https://metagenomics-workshop.readthedocs.io/en/latest/annotation/quantification.html

I performed the following steps:

1. mapped reads with bowtie2 and generated .bam file with samtools sort.
2. Removed duplicates with picard
3. Extracted gene information from prokka output files using prokkagff2gtf.sh script available from that workshop page, and obtained the .gtf file.

All seemed to go well until here.

Then I ran the following htseq command to get the read counts:

htseq-count -s no -r pos -t CDS -f bam sample2.markdup.bam sample2.map.gtf > sample2.count

The run seemed to go smoothly and produced an output file "sample2.count"

When I inspected the sample2.count, it shows gene ids in one column and second column which is supposed to show read counts, contains ALL ZERO in front of every single gene.

So, I don't know where is the problem. I would really appreciate if someone helps me in this.

Many thanks,