Add mash package

Author: fluhus

mash/mash.go

@@ -0,0 +1,51 @@
+// Package mash provides sequence MinHashing and Mash distance calculation.
+package mash
+
+import (
+	"bytes"
+	"math"
+
+	"github.com/fluhus/biostuff/sequtil/v2"
+	"github.com/fluhus/gostuff/minhash"
+	"github.com/spaolacci/murmur3"
+)
+
+// Seed is the hash seed.
+// Affects subsequent calls to [Sequences] and [Add].
+var Seed uint32 = 0
+
+// Add adds the given sequences to an existing MinHash
+// using subsequences of length k.
+// Equivalent to calling [Sequences] on the old and new sequences together.
+func Add(mh *minhash.MinHash[uint64], k int, seqs ...[]byte) {
+	h := murmur3.New64WithSeed(Seed)
+	for _, seq := range seqs {
+		for b := range sequtil.CanonicalSubsequences(bytes.ToUpper(seq), k) {
+			h.Reset()
+			h.Write(b)
+			mh.Push(h.Sum64())
+		}
+	}
+	mh.Sort()
+}
+
+// Sequences returns a single MinHash for seqs with n elements and for
+// subsequences of length k.
+func Sequences(n, k int, seqs ...[]byte) *minhash.MinHash[uint64] {
+	mh := minhash.New[uint64](n)
+	Add(mh, k, seqs...)
+	return mh
+}
+
+// Distance returns the Mash distance between two MinHash collections.
+func Distance(mh1, mh2 *minhash.MinHash[uint64], k int) float64 {
+	return FromJaccard(mh1.Jaccard(mh2), k)
+}
+
+// FromJaccard returns the Mash distance given a Jaccard similarity.
+func FromJaccard(jac float64, k int) float64 {
+	if jac == 0 {
+		return 1
+	}
+	return min(-math.Log(2*jac/(1+jac))/float64(k), 1)
+}

mash/mash_test.go

@@ -0,0 +1,50 @@
+package mash
+
+import (
+	"math"
+	"math/rand/v2"
+	"slices"
+	"testing"
+
+	"github.com/fluhus/biostuff/sequtil/v2"
+)
+
+func TestSequences(t *testing.T) {
+	seqs := [][]byte{
+		[]byte("AGATTTTTCTCCCAACGAAACTTTACAGCACGCTAGTTTACGGGCACTCC"),
+		[]byte("GCCCTATCTGGGGAGAAATCGTAGTGAGAGACCGAGGTGGCCCCACGCAC"),
+		[]byte("TACGACTGGAAGAGCCCATGCACGGATGCTGCTACTCGCATTGGTTTACG"),
+	}
+	mh1 := Sequences(5, 21, seqs...)
+	mh2 := Sequences(5, 21, seqs[0])
+	Add(mh2, 21, seqs[1])
+	Add(mh2, 21, seqs[2])
+
+	v1 := mh1.View()
+	v2 := mh2.View()
+	if !slices.Equal(v1, v2) {
+		t.Fatalf("Sequences(...) != Add(...), %v != %v", v1, v2)
+	}
+}
+
+func TestDistance(t *testing.T) {
+	const n = 1000000
+
+	// Create sequence.
+	seq1 := make([]byte, n)
+	for i := range seq1 {
+		seq1[i] = sequtil.Iton(rand.N(4))
+	}
+
+	// Create mutated sequence.
+	seq2 := slices.Clone(seq1)
+	for _, i := range rand.Perm(n)[:n/20] {
+		seq2[i] = sequtil.Iton((sequtil.Ntoi(seq2[i]) + 1 + rand.N(3)) % 4)
+	}
+
+	d := Distance(Sequences(10000, 21, seq1), Sequences(10000, 21, seq2), 21)
+	dif := math.Abs(d - 0.05)
+	if dif > 0.0001 {
+		t.Fatalf("Distance(...)=%v, want %v", d, 0.05)
+	}
+}