Hi!

A little background: I have a metagenomic sample (DNA-seq). I want to figure out the gene depth/abundance of all the genes in the sample. Since there a lot of different (and same) bacteria in the sample, I expect some of the genes to have a high abundance. Keep in mind that it is a DNA-seq, so it's not an expression analysis, but rather what genes (and their depth) the sample inhabit. Lastly, I want to group the genes in COGs while still preserving the quantitative information.

I started with an assembled contigs-file and 2 FASTQ-files (Paired-end). I then predicted genes with Prokka (bacteria only) from the contigs-file.

Next, I ran BWA where I mapped my FASTQ-files to all my predicted genes (ffn-file) and have a resulting BAM/SAM file.

Now, I need to figure out how many reads mapped to each gene. I guess I could make a small script that counts together all mapped reads, but since I'm working with paired-end data, is there some extra precautions I should take? What tools do you suggest I use? I have been told to use Featurecounts, but It doesn't seem to work, as I'm working with CDS and not Exons.

All suggestions are welcomed!