<p>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.</p>
<ol>
<li><p>Map raw reads to the contig containing your hit (original contig, pre-translation)</p>
</li>
<li><p>Normalise it with RPKM:</p>
</li>
</ol>
<p>RPKM = NumberOfReadsThatAligned / ( ContigLength/1000 * totalNumberReads/1,000,000 )</p>
<ol start="3">
<li>This will give you the relative abundance of the contig which is an APPROXIMATION of the abundance of your gene.</li>
</ol>
<p>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.</p>
<p>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.</p>
<p>Best wishes,</p>
<p>A</p>

CC BY-SA 4.0
2022-09-12T17:51:44.813Z
2022-09-12T20:41:40.907Z
2022-09-12T20:41:40.907Z
Avamys
<p>I have googled a &quot;lot&quot;, couldn't find any specific answer to the question. So, I am here seeking for your guidance. My question is similar to <a href="https://www.biostars.org/p/400176/" rel="nofollow noreferrer">this</a>. 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:</p>
<ol>
<li>Download the lasA genes from <a href="https://www.uniprot.org/uniprot/?query=gene:lasa&amp;format=fasta&amp;force=true&amp;sort=score" rel="nofollow noreferrer">uniport</a>.</li>
<li>I then made the multiple sequence alignment (lasA.msa) of the lasA gene downloaded in the previous step.</li>
<li>The protein profile was generated <code>hmmbuild res.hmm lasA.msa</code></li>
<li>Genome assembly (contigs) were done with megahit <code>megahit -1 r1.fastq.gz -2 r2.fastq.gz --min-contig-len 1000 -m 0.85 -o ASSEMBLY/ -t 20</code></li>
<li>ORFs were predicted and extracted using <code>getorf</code> from the contigs (from step 4)</li>
<li>The profile (from step 3) was used to search the ORFs using <code>hmmsearch</code></li>
<li>Several hits were obtained from <code>hmmsearch</code> indicating that the lasA gene is present in the metagenome</li>
</ol>
<p>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.</p>
<p>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?</p>
<p>Thanks and regards</p>

<genome><assembly><metagenome><bowtie2><hmmer>
How to quantifiy of specific genes from shotgun metagenome?
# How to quantifiy of specific genes from shotgun metagenome?

(Asked by: username_0 on 2021-07-18)

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](https://www.biostars.org/p/400176/). 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](https://www.uniprot.org/uniprot/?query=gene:lasa&format=fasta&force=true&sort=score).
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

**Tags:** genome, assembly, metagenome, bowtie2, hmmer

---

## 1 answers

### Answer by username_1 on 2022-09-12. Score: 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


---

