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:

  1. Download the lasA genes from uniport.
  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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.