How to quantifiy of specific genes from shotgun metagenome?

Question

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

Answer 1

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:

RPKM = NumberOfReadsThatAligned / ( ContigLength/1000 * totalNumberReads/1,000,000 )

3. This will give you the relative abundance of the contig which is an APPROXIMATION of the abundance of your gene.

Alternatively, to be more precise you can FragGeneScan your raw reads, extract the sequence of your putative protein and map translated raw reads to it.

For your other question: Yes, HMMs are more sensitive than blast or diamond since its output is probability score that these protein contain motifs extracted from your alignments (given that there is one). For that reason you need to look at the bitscore and E-value. If the bitscore is low and/or E-value is high, you have a high chance of false-positive.

Best wishes,

A