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:
- Download the lasA genes from uniport.
- I then made the multiple sequence alignment (lasA.msa) of the lasA gene downloaded in the previous step.
- The protein profile was generated
hmmbuild res.hmm lasA.msa
- 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
- ORFs were predicted and extracted using
getorf
from the contigs (from step 4) - The profile (from step 3) was used to search the ORFs using
hmmsearch
- 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