I am using the following metagenomic workshop tutorial to analyse my own metagenomic data.
I performed the following steps:
- mapped reads with bowtie2 and generated .bam file with samtools sort.
- Removed duplicates with picard
- Extracted gene information from prokka output files using
prokkagff2gtf.shscript 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.