# How to quantifiy of specific genes from shotgun metagenome?

I have googled a "lot", couldn't find any specific answer to the question. So, I am here seeking for your guidance. My question is similar to this. I have several metagenome (n=30). But for simplicity’s sake, let's assume only one metagenome. I wanted to quantify lasA genes in the metagenome. So far, what I have done (not sure if it's the right approach) is:

2. I then made the multiple sequence alignment (lasA.msa) of the lasA gene downloaded in the previous step.
3. The protein profile was generated hmmbuild res.hmm lasA.msa
4. Genome assembly (contigs) were done with megahit megahit -1 r1.fastq.gz -2 r2.fastq.gz --min-contig-len 1000 -m 0.85 -o ASSEMBLY/ -t 20
5. ORFs were predicted and extracted using getorf from the contigs (from step 4)
6. The profile (from step 3) was used to search the ORFs using hmmsearch
7. Several hits were obtained from hmmsearch indicating that the lasA gene is present in the metagenome

Now, my question is, how do I count the number of reads affiliated to the lasA gene? I feel that I need to map the reads to the assembly, but I don't know how to go about that.

NOTE: I have used MG-RAST, Humann3, and eggNOG-mapper, but none of them showed the presence of lasA gene, however the above-mentioned hmmsearch was able to find the gene, what could be the reason? Problem with my approach? Or hmmsearch using the protein profile is more sensitive?

Thanks and regards

The way I would do it is to make an assumption that the reads mapping to your contig == (abundance of the genes in the contig)/contig length.

1. Map raw reads to the contig containing your hit (original contig, pre-translation)

2. Normalise it with RPKM: