using sourmash to select best genome for mapping

[ctb]

This came up on microbial bioinformatics slack, thought I'd share - topic was selecting from many viral genomes.

it should be possible to use large-scale ANI-style analyses to select the closest genome for mapping. we’ve been doing this with sourmash and genome-grist for metagenomes, and I know tools like ganon and I think kmcp can do the same thing.
with sourmash I would say the first thing to try is:

• sketch all your genome references with sourmash sketch dna -p scaled=100 *.fna
• do the same with your metagenome(s)/shotgun reads
• run sourmash prefetch <metagenome>.sig genome*.fna.sig --threshold-bp=0 -o matches.csv
• sort matches.csv on f_match_query and pick the highest value (this is "k-mer detection" per use 'detection' terminology for fraction-of-genome-kmers-found #2170) and use that match_name as the reference genome.

happy to help troubleshoot here or on sourmash issue tracker if there is interest.

A more fun and “sophisticated” approach that could go horribly awry is to use sourmash gather after the prefetch to develop a minimum metagenome cover, but that’s only for people that are ok with expending some of their time and energy on a potential wild goose chase (which I am happy to support, but, you know, still)

[ctb]

Added to #2184 in 69dfbc0:

Can I use sourmash to determine the best reference genome for mapping my reads?

Yes! (And see the FAQ above,
How do k-mer analyses compare with read mapping?

If you're interested in picking a single best reference genome (from a
large database) for read mapping, you can do the following:

• sketch all your reference genomes with sourmash sketch dna, and/or download
  one of our prepared databases. You can use default parameters for sourmash sketch.
• sketch your read collection with sourmash sketch dna.
• run sourmash prefetch reads.sig <reference genome sketch collection(s)> -o matches.csv
• sort matches.csv on the f_match_query column and pick the highest value -
  this is the k-mer detection - and pick the match name from the match_name column;
• use that reference genome as your mapping target.

If you want to map a metagenome to multiple references, consider
using sourmash gather and/or
the genome-grist workflow.

(This is also known as "read recruitment.")