<p>Thats an interesting question. I do this stuff all the time and don't have a fixed cut-off. What I am certain about is 100 is marginal evidence of a species. So I personally would be reluctant to trust this - albeit my stuff is not gut microbes. I would consider approximately one thousand as the threshold of a species designation (see below), 500 at a minimum. The &quot;copy numbers&quot;, i.e. number of hits per genome, can be extraordinary. Thus species with 100 usually drops off the map. The rationale is simple, there are upwards of 50 million records per data set.</p>
<p>I fully understand the basis of the question. Metagenomes once the paired end are assembled (not mine) can result in log10 reduction in total number of reads and this can cause problems of sample size (just check the number of reads before and after assembly). Under those circumstances 100 might be okay because the large reduction in sample would be reflected in a reduced number of hits.</p>
<p>I cannot recommend 20 as evidence of a species. Artefacts can happen in metagenomics. The genomes of the Kraken2 database are never a perfect match to the in situ genetic material, and Kraken2 accommodates that, however errors, i.e. non-specific binding, are inevitable. There are millions of reads and non-specific binding absolutely happen. To prove it take a 20 metagenome reads and blast them - you'll find non-specific hits.</p>
<p>Twenty random hits for a species is entirely feasible in my opinion specifically for a metagenomics data set (it's a special case). I'm not aware of anyone who has empirically investigated this, but it would be a difficult investigation.</p>
<hr />
<p>Just to be clear this question is specifically about metagenomics. Investigators use Kraken2 for all sorts of things because it's an excellent algorithm. The above doesn't apply to the &quot;all sorts of things&quot; stuff.</p>
<hr />
<p>Using Brocken as @gringer is preferred but for presence/absence of species (not abundance estimation) Kraken2 alone is fine.</p>

CC BY-SA 4.0
2024-06-21T11:57:15.653Z
Thank you for your detailed answer! I agree that 100 reads for species is not a good cutoff. I'm using an assembled isolate of ours added to the custom db of Kraken. I want to compare if this isolate is in patient datasets rather than healthy datasets. I can compare abundance differences without issues but saying if the specific strain isolate we have is in the samples 1/22/65... is getting a bit tricky. Maybe I can say something like ~10.000 reads are mapped to E.coli which proves species are there and 100 of them mapped to our isolate X which is a good indicator that the strain is there?
ahmet
CC BY-SA 4.0
2024-06-21T13:49:07.783Z
@ahmet upvotes and "accepts" (green tick) are strongly preferred to notes of thanks. What you are describing is a separate question and the site format is one question per post.
M__
CC BY-SA 4.0
2024-06-21T14:23:07.817Z
ah, sorry I thought it was the same root question but a bit elaborated. Unfortunately, the site does not allow me to upvote yet because of reputation limits. Thanks anyway!
ahmet
CC BY-SA 4.0
2024-06-21T14:42:40.787Z
Huge thanks @ahmed. If you submit this as a separate question I'll forward an answer and it will provide you with the threshold to upvote in the future.
M__
CC BY-SA 4.0
2024-06-24T08:16:07.853Z
What do you mean with "non-specific binding" in this context?
Pallie
CC BY-SA 4.0
2024-06-24T16:23:16.963Z
@Pallie false positives
M__
CC BY-SA 4.0
2024-06-28T10:44:06.587Z
@M__In the sense of the sequencing read produced by the instrument or in the sense of Kraken2 misclassification?
Pallie
CC BY-SA 4.0
2024-06-21T04:47:06.950Z
2024-06-21T13:52:25.930Z
2024-06-21T13:52:25.930Z
M__
<p>Given that you're already using kraken2, I recommend bracken:</p>
<p><a href="https://github.com/jenniferlu717/Bracken" rel="nofollow noreferrer">https://github.com/jenniferlu717/Bracken</a></p>
<blockquote>
<p>Bracken is a companion program to Kraken 1, KrakenUniq, or Kraken 2 While Kraken classifies reads to multiple levels in the taxonomic tree, Bracken allows estimation of abundance at a single level using those classifications (e.g. Bracken can estimate abundance of species within a sample).</p>
</blockquote>
<p>One thing to bear in mind with kmer-based approaches for metagenomics, especially when short reads are involved, is that there is a substantial false positive hit rate. In other words, don't trust the result of a mapping to be 100% accurate. This is even more important when considering that horizontal gene transfer is fairly common in bacteria; some people may not believe a bacteria is definitely present until they have seen a fully-assembled genome.</p>
<p>As an aside, because sequencing is an incomplete, [hopefully] random sampling, it's impossible to say with full confidence that a particular bacterium is <em>not</em> present in a sample. It may just be that the sampling method (i.e. sample prep) didn't select the bacterium for sequencing. Given that, it's best to think about limits of detection.</p>
<p>My general rule of thumb is to sample about ten times the minimum amount required for detection. Detecting down to 1% presence in the sample with perfect read classification and no false positive/negative assignment, I would want to see at least 1,000 reads. For 0.1%, I'd want to see at least 10,000 reads. To account for false positives and negatives, I might increase that by another factor of 10 to be safe (e.g. 100,000 reads for a 0.1% detection limit). In these cases, the statements I would make would be of the form, &quot;Bacterium X is unlikely to be present in the sample at a level greater than y%.&quot;</p>
<p>These answers are <em>always</em> going to be project-specific. It is up to the analysts, researchers, and biologists to decide where to put the thresholds. There will never be a universal threshold that works in every situation.</p>
<p>And finally, never put complete trust in a single technology or method. Try to validate results with different technologies to reduce the risk of systematic error.</p>

CC BY-SA 4.0
2024-06-21T12:03:53.763Z
Thank you for your detailed answer! Unfortunately Bracken gives errors when i try to compile my custom db with few added isolates of ours. And the developers don't get back to the open issues on github nowadays. But I'll definitely try again. I guess in my case (please see my comment on M__'s answer) I can investigate the read counts of specific cases between patient and healthy datasets with different conditions and set some sort of confidence thresholds using 16s and short read map percentages. And hopefully see some trends fitting our hypothesis.
ahmet
CC BY-SA 4.0
2024-06-21T22:25:49.960Z
If you've got errors with running Bracken, feel free to ask a separate question about that.
gringer
CC BY-SA 4.0
2024-06-21T05:40:47.523Z
2024-06-21T05:40:47.523Z
gringer
<p>I have ~100 paired end short read human gut metagenome samples that I classified using Kraken2. Now I want to know if a specific bacterium is in any of those samples. As far as I've searched people use mapped read thresholds ranging from 20-100 reads or percentage of the reads mapped to a specific bacterium to the entire mapped bacterial reads around 0.001 to 0.01. Or a mix of coverage and depth. But I didn't see any generic method to find how many reads should be mapped to a reference genome to say that some bacterium is in a metagenome.</p>
<p>Is there a best practice sort of way to determine if a bacterium is in some metagenomic sample?</p>

CC BY-SA 4.0
2024-06-21T14:51:31.373Z
Until I have the time to write a full answer, I suggest you look into YACHT: https://pubmed.ncbi.nlm.nih.gov/37131762/
acvill
CC BY-SA 4.0
2024-06-21T15:03:41.220Z
ah, I didn't check the preprints! That might be exactly what I'm looking for, thanks!
ahmet
<r><illumina><metagenome><bacteria><taxonomy>
short Read/percentage threshold for bacterium presence in metagenome
# short Read/percentage threshold for bacterium presence in metagenome

(Asked by: username_0 on 2024-06-20)

I have ~100 paired end short read human gut metagenome samples that I classified using Kraken2. Now I want to know if a specific bacterium is in any of those samples. As far as I've searched people use mapped read thresholds ranging from 20-100 reads or percentage of the reads mapped to a specific bacterium to the entire mapped bacterial reads around 0.001 to 0.01. Or a mix of coverage and depth. But I didn't see any generic method to find how many reads should be mapped to a reference genome to say that some bacterium is in a metagenome.

Is there a best practice sort of way to determine if a bacterium is in some metagenomic sample?

**Tags:** r, illumina, metagenome, bacteria, taxonomy
> **Comments:**
> **username_1**: Until I have the time to write a full answer, I suggest you look into YACHT: https://pubmed.ncbi.nlm.nih.gov/37131762/  
> **username_0**: ah, I didn't check the preprints! That might be exactly what I'm looking for, thanks!  

---

## 2 answers

### Answer by username_2 on 2024-06-21. Score: 2 | Accepted answer

Thats an interesting question. I do this stuff all the time and don't have a fixed cut-off. What I am certain about is 100 is marginal evidence of a species. So I personally would be reluctant to trust this - albeit my stuff is not gut microbes. I would consider approximately one thousand as the threshold of a species designation (see below), 500 at a minimum. The "copy numbers", i.e. number of hits per genome, can be extraordinary. Thus species with 100 usually drops off the map. The rationale is simple, there are upwards of 50 million records per data set.

I fully understand the basis of the question. Metagenomes once the paired end are assembled (not mine) can result in log10 reduction in total number of reads and this can cause problems of sample size (just check the number of reads before and after assembly). Under those circumstances 100 might be okay because the large reduction in sample would be reflected in a reduced number of hits.

I cannot recommend 20 as evidence of a species. Artefacts can happen in metagenomics. The genomes of the Kraken2 database are never a perfect match to the in situ genetic material, and Kraken2 accommodates that, however errors, i.e. non-specific binding, are inevitable. There are millions of reads and non-specific binding absolutely happen. To prove it take a 20 metagenome reads and blast them - you'll find non-specific hits.

Twenty random hits for a species is entirely feasible in my opinion specifically for a metagenomics data set (it's a special case). I'm not aware of anyone who has empirically investigated this, but it would be a difficult investigation.

---

Just to be clear this question is specifically about metagenomics. Investigators use Kraken2 for all sorts of things because it's an excellent algorithm. The above doesn't apply to the "all sorts of things" stuff.

---

Using Brocken as @gringer is preferred but for presence/absence of species (not abundance estimation) Kraken2 alone is fine.

> **Comments:**
> **username_0**: Thank you for your detailed answer! I agree that 100 reads for species is not a good cutoff. I'm using an assembled isolate of ours added to the custom db of Kraken. I want to compare if this isolate is in patient datasets rather than healthy datasets. I can compare abundance differences without issues but saying if the specific strain isolate we have is in the samples 1/22/65... is getting a bit tricky. Maybe I can say something like ~10.000 reads are mapped to E.coli which proves species are there and 100 of them mapped to our isolate X which is a good indicator that the strain is there?  
> **username_2**: @ahmet upvotes and "accepts" (green tick) are strongly preferred to notes of thanks. What you are describing is a separate question and the site format is one question per post.  
> **username_0**: ah, sorry I thought it was the same root question but a bit elaborated. Unfortunately, the site does not allow me to upvote yet because of reputation limits. Thanks anyway!  
> **username_2**: Huge thanks @ahmed. If you submit this as a separate question I'll forward an answer and it will provide you with the threshold to upvote in the future.  
> **username_3**: What do you mean with "non-specific binding" in this context?  
> **username_2**: @Pallie false positives  
> **username_3**: @M__In the sense of the sequencing read produced by the instrument or in the sense of Kraken2 misclassification?  

---

### Answer by username_4 on 2024-06-21. Score: 2

Given that you're already using kraken2, I recommend bracken:

<https://github.com/jenniferlu717/Bracken>

> Bracken is a companion program to Kraken 1, KrakenUniq, or Kraken 2 While Kraken classifies reads to multiple levels in the taxonomic tree, Bracken allows estimation of abundance at a single level using those classifications (e.g. Bracken can estimate abundance of species within a sample).

One thing to bear in mind with kmer-based approaches for metagenomics, especially when short reads are involved, is that there is a substantial false positive hit rate. In other words, don't trust the result of a mapping to be 100% accurate. This is even more important when considering that horizontal gene transfer is fairly common in bacteria; some people may not believe a bacteria is definitely present until they have seen a fully-assembled genome.

As an aside, because sequencing is an incomplete, [hopefully] random sampling, it's impossible to say with full confidence that a particular bacterium is *not* present in a sample. It may just be that the sampling method (i.e. sample prep) didn't select the bacterium for sequencing. Given that, it's best to think about limits of detection.

My general rule of thumb is to sample about ten times the minimum amount required for detection. Detecting down to 1% presence in the sample with perfect read classification and no false positive/negative assignment, I would want to see at least 1,000 reads. For 0.1%, I'd want to see at least 10,000 reads. To account for false positives and negatives, I might increase that by another factor of 10 to be safe (e.g. 100,000 reads for a 0.1% detection limit). In these cases, the statements I would make would be of the form, "Bacterium X is unlikely to be present in the sample at a level greater than y%."

These answers are *always* going to be project-specific. It is up to the analysts, researchers, and biologists to decide where to put the thresholds. There will never be a universal threshold that works in every situation.

And finally, never put complete trust in a single technology or method. Try to validate results with different technologies to reduce the risk of systematic error.

> **Comments:**
> **username_0**: Thank you for your detailed answer! Unfortunately Bracken gives errors when i try to compile my custom db with few added isolates of ours. And the developers don't get back to the open issues on github nowadays. But I'll definitely try again. I guess in my case (please see my comment on M__'s answer) I can investigate the read counts of specific cases between patient and healthy datasets with different conditions and set some sort of confidence thresholds using 16s and short read map percentages. And hopefully see some trends fitting our hypothesis.  
> **username_4**: If you've got errors with running Bracken, feel free to ask a separate question about that.  

---

