Minimap2 detects too many 16S sequences in metagenome

Question

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.

Answer 1

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.

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 (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 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:

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.