MinHash Sketch is a method of rapidly comparing large strings or sets. In genomics, you can use it like this:
1) Gather all the kmers in a genome.
2) Apply a hash function to them.
3) Keep the 10000 smallest hashcodes and call this set a “sketch”.
If you do this for multiple genomes, you can calculate how similar two genomes are much faster than via alignment; and specifically, the greatest advantage is that the speed of sketch comparison is unrelated to genome size. For example, if you assemble something, you can sketch it and compare it to sketches of everything in nt or RefSeq in a couple seconds. If the top hit shows that your new sketch shared 90% of its hash codes with E.coli, that means it has roughly 90% kmer identity to E.coli, and therefore, it is E.coli.
The BBMap package has an extremely efficient MinHash Sketch implementation, now accessible through 4 programs. The latest, just released in BBMap 36.92, is SendSketch:
sendsketch.sh in=contigs.fa
This will make a sketch of your assembly, open an HTTP connection to JGI’s taxonomy server, and send the sketch. The server will compare it to sketches of all of nt, and return the top hits. This is kind of convenient because the sketches sit around in memory all the time rather than needing to be loaded, so you get the answer in a couple seconds. You can also add the flag “refseq” to compare to RefSeq instead of nt. Additionally, SendSketch works on raw reads (including Illumina, PacBio, and Nanopore). With Illumina the files can be huge and you don’t need all of the reads so I normally run it like this:
sendsketch.sh in=reads.fq reads=400k
…to limit it to the first 400,000 reads. Raw PacBio has a high error rate so I typically increase the default sketch size to increase sensitivity:
sendsketch.sh in=PacBio.fq size=100000
You can also run a local server with “taxserver.sh” if you want to use a different database for SendSketch.
The other tools are for processing sketches locally. To sketch and query refseq, for example:
sketch.sh in=refseq.fa.gz out=refseq#.sketch files=31 mode=sequence
That will produce one sketch per sequence, and split them among 31 output files (for rapid loading). Alternatively, and preferably, you can produce one sketch per taxID. This process is a bit more involved, but the full details are given in the BBMap package in the /bbmap/pipelines/ directory, in these example shellscripts:
fetchtaxonomy.sh
fetchnt.sh
fetchRefSeq.sh
Now that you have created some sketches, you can query them like this:
comparesketch.sh in=contigs.fa refseq*.sketch
“in=” can be a comma-delimited list of either fasta or sketch files. Also, you can turn an input fasta into one sketch per sequence rather than just a single sketch with the “mode=sequence” flag. If you want to process with CompareSketch using taxonomic information, add “tree=taxtree.gz” as well (which was created by the “fetchtaxonomy.sh” script). SendSketch automatically uses taxonomy on the server-side.
So, that’s how you run Sketch. Now, what does the output mean? Here’s an example, using E.coli:
The complete description of columns is in /bbmap/docs/guides/BBSketchGuide.txt, but briefly, “matches” is the number of matching kmers between the query and reference sketch. KID (or kmer identity) is simply the percent of the kmers that matched, while WKID is the weighted kmer identity, which has been adjusted to compensate for genome size. The ANI, completeness, and contamination are estimates, comparing the query to that reference.
Titus Brown suggested that I cite prior art, since this is not a novel idea, so:
The MinHash concept, which is about 20 years old: en.wikipedia.org/wiki/MinHash
Mash: github.com/marbl/Mash
SourMash: github.com/dib-lab/sourmash
Read more here: Source link
