Answers.Body: <p>16S is a very conserved sequence, which is why it's used for targeted phylogenetic analysis; it makes it easy to amplify and analyse. Unfortunately that conservation is an issue with minimap2, which is built around the idea that matching scattered subsequences within a sequence is good enough for identifying matches to the entire sequence. With a highly-conserved sequence, there's a higher chance that subsequences will be identical, meaning that diverged sequences will appear to be the identical from the perspective of the mapper.</p>
  <p>Minimap2 can be tweaked to make it less sensitive and more specific, but you'd be better off using a different mapping tool, especially one designed for phylogenetic analysis and conserved sequences. I recommend <a href="https://github.com/DerrickWood/kraken2/wiki/Manual#installation" rel="nofollow noreferrer">kraken2</a>.</p>
  <p>pre-built 16S indexes for kraken2 can be found here:</p>
  <p><a href="https://benlangmead.github.io/aws-indexes/k2" rel="nofollow noreferrer">https://benlangmead.github.io/aws-indexes/k2</a></p>
  <p>If you really want to map specifically to your own generated metagenome, then you could try <a href="https://gitlab.com/mcfrith/last" rel="nofollow noreferrer">LAST</a> (with results filtered through <code>last-map-probs</code>), which I have found to be a little bit better in dealing with specificity issues.</p>
  <p>Aside: I notice that you've stated &quot;bacterial sequences from a <em>metagenome</em>&quot;. If these are, in fact, metagenomic reads (and not targeted 16S reads), then you would be better off using one of the larger metagenome bacterial databases  (e.g. plusPF or plusPF-8, depending on computer memory) instead of concentrating only on 16S-matching reads. This will provide much better resolution; 16S only really works down to the family level.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-04-20T01:51:03.400Z
Answers.LastActivityDate: 2022-04-20T01:51:03.400Z
Answers.UserDisplayName: gringer
Answers.Body: <p>You are taking reads as the reference. Given the similarity of 16S, you will find hits for most 16S in other species. You have to filter out poor hits. This is the faster way if you don't have many contigs.</p>
  <p>Alternatively, switch the order of reference and reads. But in this case, you need to use a high k-mer occurrence cut off as most of your k-mers will occur &gt;10,000 times. Something like:</p>
  <pre class="lang-sh prettyprint-override"><code>minimap2 -c -f30000 -N100 ref-16S.fa reads.fq
  </code></pre>
  <p>This will be fairly slow as there are too many hits in the database but it will be easier to process the output.</p>
  <p>You can try both and see what works for you.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-04-22T21:50:52.683Z
Answers.LastActivityDate: 2022-04-22T21:50:52.683Z
Answers.UserDisplayName: user172818
Body: <p>I'd like to extract bacterial sequences from a metagenome where ~99% of contigs are from a host insect. My current protocol minimaps every 16S sequences present in the sample, then blasts all those sequences to create a library of full genomes that I minimap against. I'm using the full metagenome as my reference file. I wasn't sure if that makes a difference, but figured I'd mention it. I ran the following command to retrieve 16S sequences that are present in the sample:</p>
  <pre><code>minimap2 All_reads.fastq 16S_reference_library.fasta &gt; 16S_map.paf
  </code></pre>
  <p>Minimap2 retrieves every entry in 16S_reference_library.fasta. It's highly improbable some of those sequences are present, or that 25,000 species are present at all. Are there any parameters I could adjust in Minimap2 that might help with this?</p>
  <p>Alternatively, I'm trying to implement Krakenuniq and Centrifuge, but I'm concerned that, because I don't have a reference genome for the host insect, they'll have a very high false positive rate.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-04-19T23:12:55.333Z
Comments.Text: Do you mean to be using your reads (`fastq` file) as the reference rather than the query? Normal syntax is `minimap2 reference.fasta query.fastq`. What do you mean by "retrieves"? How good are the alignments, are they `*`? Can you show us a snippet of the PAF?
Comments.UserDisplayName: Maximilian Press
Tags: <metagenome><minimap2>
Title: Minimap2 detects too many 16S sequences in metagenome
ThreadText: # Minimap2 detects too many 16S sequences in metagenome
  
  (Asked by: username_0 on 2022-04-19)
  
  I'd like to extract bacterial sequences from a metagenome where ~99% of contigs are from a host insect. My current protocol minimaps every 16S sequences present in the sample, then blasts all those sequences to create a library of full genomes that I minimap against. I'm using the full metagenome as my reference file. I wasn't sure if that makes a difference, but figured I'd mention it. I ran the following command to retrieve 16S sequences that are present in the sample:
  
  ```
  minimap2 All_reads.fastq 16S_reference_library.fasta > 16S_map.paf
  
  ```
  
  Minimap2 retrieves every entry in 16S\_reference\_library.fasta. It's highly improbable some of those sequences are present, or that 25,000 species are present at all. Are there any parameters I could adjust in Minimap2 that might help with this?
  
  Alternatively, I'm trying to implement Krakenuniq and Centrifuge, but I'm concerned that, because I don't have a reference genome for the host insect, they'll have a very high false positive rate.
  
  **Tags:** metagenome, minimap2
  
  > **Comments:**
  > **username_1**: Do you mean to be using your reads (`fastq` file) as the reference rather than the query? Normal syntax is `minimap2 reference.fasta query.fastq`. What do you mean by "retrieves"? How good are the alignments, are they `*`? Can you show us a snippet of the PAF?  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2022-04-20. Score: 3
  
  16S is a very conserved sequence, which is why it's used for targeted phylogenetic analysis; it makes it easy to amplify and analyse. Unfortunately that conservation is an issue with minimap2, which is built around the idea that matching scattered subsequences within a sequence is good enough for identifying matches to the entire sequence. With a highly-conserved sequence, there's a higher chance that subsequences will be identical, meaning that diverged sequences will appear to be the identical from the perspective of the mapper.
  
  Minimap2 can be tweaked to make it less sensitive and more specific, but you'd be better off using a different mapping tool, especially one designed for phylogenetic analysis and conserved sequences. I recommend [kraken2](https://github.com/DerrickWood/kraken2/wiki/Manual#installation).
  
  pre-built 16S indexes for kraken2 can be found here:
  
  <https://benlangmead.github.io/aws-indexes/k2>
  
  If you really want to map specifically to your own generated metagenome, then you could try [LAST](https://gitlab.com/mcfrith/last) (with results filtered through `last-map-probs`), which I have found to be a little bit better in dealing with specificity issues.
  
  Aside: I notice that you've stated "bacterial sequences from a *metagenome*". If these are, in fact, metagenomic reads (and not targeted 16S reads), then you would be better off using one of the larger metagenome bacterial databases (e.g. plusPF or plusPF-8, depending on computer memory) instead of concentrating only on 16S-matching reads. This will provide much better resolution; 16S only really works down to the family level.
  
  
  ---
  
  ### Answer by username_3 on 2022-04-22. Score: 2
  
  You are taking reads as the reference. Given the similarity of 16S, you will find hits for most 16S in other species. You have to filter out poor hits. This is the faster way if you don't have many contigs.
  
  Alternatively, switch the order of reference and reads. But in this case, you need to use a high k-mer occurrence cut off as most of your k-mers will occur >10,000 times. Something like:
  
  ```sh
  minimap2 -c -f30000 -N100 ref-16S.fa reads.fq
  
  ```
  
  This will be fairly slow as there are too many hits in the database but it will be easier to process the output.
  
  You can try both and see what works for you.
  
  
  ---
  
  