diff --git a/docmd/implementation/obicompactvec.md b/docmd/implementation/obicompactvec.md index 301b021..34c2024 100644 --- a/docmd/implementation/obicompactvec.md +++ b/docmd/implementation/obicompactvec.md @@ -347,11 +347,24 @@ Provided finalisations: | `relfreq_euclidean_dist_matrix()` | `√partial_relfreq_euclidean[i,j]` | | `hellinger_dist_matrix()` | `√partial_hellinger[i,j] / √2` | | `hellinger_euclidean_dist_matrix()` | `√partial_hellinger[i,j]` | +| `threshold_mash_dist_matrix(k, t)` | Mash distance, derived from `threshold_jaccard_dist_matrix(t)` — no separate partial | ### BitPartials Required: `partial_jaccard() -> (Array2, Array2)`, `partial_hamming() -> Array2`. Both additive across layers and partitions. +Provided finalisations also include `jaccard_dist_matrix()`, `hamming_dist_matrix()`, and `mash_dist_matrix(k)`. + +### Mash distance + +`mash_dist_matrix`/`threshold_mash_dist_matrix` add no new additive primitive: both are a pointwise transform of the existing Jaccard distance matrix, per the Mash mutation-rate estimator [@Mash-distances-doc; @Fan2015-mash-formula]: + +``` +D = -1/k · ln(2J / (1+J)), J = 1 - d_jaccard +``` + +`J ≤ 0` (i.e. `d_jaccard ≥ 1`, no shared k-mers) maps to `D = 1` (maximal distance) rather than the `ln` singularity at `J = 0`. + --- ## Temp-file-backed types diff --git a/docmd/index.md b/docmd/index.md index 9f6d650..812cc86 100644 --- a/docmd/index.md +++ b/docmd/index.md @@ -13,7 +13,7 @@ | `query` | Query an index with sequences and annotate matches | | `dump` | Dump all indexed k-mers as CSV (kmer + per-genome counts or presence); supports the shared [kmer filtering](implementation/filtering.md) system; `--head N` limits output to the first N k-mers | | `annotate` | Add or update genome metadata from a CSV file; or dump metadata as CSV | -| `distance` | Compute pairwise distance matrix between genomes; optionally build NJ/UPGMA trees; `--presence-threshold N` sets the minimum count to consider a k-mer present when computing Jaccard on count indexes (default 1) | +| `distance` | Compute pairwise distance matrix between genomes (`--metric jaccard\|mash\|hamming\|bray-curtis\|relfreq-bray-curtis\|euclidean\|relfreq-euclidean\|hellinger\|hellinger-euclidean`); optionally build NJ/UPGMA trees; `--presence-threshold N` sets the minimum count to consider a k-mer present when computing Jaccard/Mash on count indexes (default 1) | | `unitig` | Build a global de Bruijn graph across all partitions and enumerate its unitigs as FASTA; supports the shared [kmer filtering](implementation/filtering.md) system | | `select` | Project and/or aggregate genome columns into a new or in-place index; the column-axis counterpart of `filter` (see [select](implementation/select.md)) | | `estimate` | Estimate approximate-index parameters (z, evidence bits, FP rates) before indexing | diff --git a/docmd/references.bib b/docmd/references.bib index a195b12..fca4c8b 100644 --- a/docmd/references.bib +++ b/docmd/references.bib @@ -241,3 +241,21 @@ volume = 33, year = 2017, bdsk-url-1 = {http://dx.doi.org/10.1093/bioinformatics/btw832}} + +@misc{Mash-distances-doc, + author = {{Marbl Lab}}, + howpublished = {Mash documentation}, + title = {Mash Distance}, + url = {https://mash.readthedocs.io/en/latest/distances.html}, + urldate = {2026-07-09}, + year = 2026} + +@article{Fan2015-mash-formula, + author = {Fan, Huan and Ives, Anthony R and Surget-Groba, Yann and Cannon, Charles H}, + doi = {10.1186/s12864-015-1647-5}, + journal = {BMC Genomics}, + number = 1, + title = {An assembly and alignment-free method of phylogeny reconstruction from next-generation sequencing data}, + url = {https://doi.org/10.1186/s12864-015-1647-5}, + volume = 16, + year = 2015} diff --git a/src/obicompactvec/src/traits.rs b/src/obicompactvec/src/traits.rs index cc52bc1..f3c0440 100644 --- a/src/obicompactvec/src/traits.rs +++ b/src/obicompactvec/src/traits.rs @@ -1,5 +1,16 @@ use ndarray::{Array1, Array2}; +/// Convert a Jaccard distance matrix (`1 - J`) into a Mash distance matrix, per +/// https://mash.readthedocs.io/en/latest/distances.html: +/// `D = -1/k * ln(2J / (1+J))`. +fn jaccard_to_mash(d_jaccard: &Array2, k: usize) -> Array2 { + d_jaccard.mapv(|d| { + let j = 1.0 - d; + if j <= 0.0 { 1.0 } + else { -1.0 / k as f64 * (2.0 * j / (1.0 + j)).ln() } + }) +} + // ── Column-level weight statistic — total count or presence count per column. /// Additive across layers and partitions; used as denominator in normalised distances. /// @@ -74,6 +85,12 @@ pub trait CountPartials: ColumnWeights { m } + /// Mash distance (https://mash.readthedocs.io/en/latest/distances.html), derived + /// from the presence-threshold Jaccard distance. + fn threshold_mash_dist_matrix(&self, k: usize, threshold: u32) -> Array2 { + jaccard_to_mash(&self.threshold_jaccard_dist_matrix(threshold), k) + } + fn relfreq_bray_dist_matrix(&self) -> Array2 { let global = self.col_weights(); let mut m = self.partial_relfreq_bray(&global).mapv(|v| 1.0 - v); @@ -126,6 +143,12 @@ pub trait BitPartials: ColumnWeights { m } + /// Mash distance (https://mash.readthedocs.io/en/latest/distances.html), derived + /// from the Jaccard distance. + fn mash_dist_matrix(&self, k: usize) -> Array2 { + jaccard_to_mash(&self.jaccard_dist_matrix(), k) + } + fn hamming_dist_matrix(&self) -> Array2 { self.partial_hamming() } diff --git a/src/obikindex/src/distance.rs b/src/obikindex/src/distance.rs index 867cff0..f5d63be 100644 --- a/src/obikindex/src/distance.rs +++ b/src/obikindex/src/distance.rs @@ -14,6 +14,8 @@ pub enum DistanceMetric { Jaccard, /// Hamming distance (number of differing kmer positions) on presence/absence data. Hamming, + /// Mash distance on presence/absence data (Jaccard-derived mutation-rate estimate). + Mash, /// Bray-Curtis dissimilarity on raw counts. BrayCurtis, /// Bray-Curtis dissimilarity normalised by per-genome total counts. @@ -84,6 +86,7 @@ impl KmerIndex { DistanceMetric::Hellinger => CountPartials::hellinger_dist_matrix(&global), DistanceMetric::HellingerEuclidean => CountPartials::hellinger_euclidean_dist_matrix(&global), DistanceMetric::Jaccard => CountPartials::threshold_jaccard_dist_matrix(&global, presence_threshold), + DistanceMetric::Mash => CountPartials::threshold_mash_dist_matrix(&global, self.kmer_size(), presence_threshold), DistanceMetric::Hamming => { return Err(OKIError::InvalidInput( "Hamming is only available for presence/absence indexes".into(), @@ -108,6 +111,7 @@ impl KmerIndex { let matrix = match metric { DistanceMetric::Jaccard => BitPartials::jaccard_dist_matrix(&global), + DistanceMetric::Mash => BitPartials::mash_dist_matrix(&global, self.kmer_size()), DistanceMetric::Hamming => { BitPartials::hamming_dist_matrix(&global).mapv(|v| v as f64) } diff --git a/src/obikmer/src/cmd/distance.rs b/src/obikmer/src/cmd/distance.rs index f66d78f..999beba 100644 --- a/src/obikmer/src/cmd/distance.rs +++ b/src/obikmer/src/cmd/distance.rs @@ -10,6 +10,7 @@ use tracing::info; #[derive(clap::ValueEnum, Clone, Copy, Debug)] pub enum MetricArg { Jaccard, + Mash, Hamming, BrayCurtis, #[value(name = "relfreq-bray-curtis")] @@ -26,6 +27,7 @@ impl From for DistanceMetric { fn from(m: MetricArg) -> Self { match m { MetricArg::Jaccard => DistanceMetric::Jaccard, + MetricArg::Mash => DistanceMetric::Mash, MetricArg::Hamming => DistanceMetric::Hamming, MetricArg::BrayCurtis => DistanceMetric::BrayCurtis, MetricArg::RelfreqBrayCurtis => DistanceMetric::RelfreqBrayCurtis, diff --git a/src/obiskbuilder/src/rolling_stat.rs b/src/obiskbuilder/src/rolling_stat.rs index 09b3f8d..e8f1bc1 100644 --- a/src/obiskbuilder/src/rolling_stat.rs +++ b/src/obiskbuilder/src/rolling_stat.rs @@ -196,13 +196,6 @@ impl RollingStat { .map(|raw| Minimizer::from_raw_unchecked(raw << (64 - self.m * 2))) } - pub fn entropy(&self, order: usize) -> Option { - if !self.ready() { - return None; - } - Some(self.entropy.entropy(order)) - } - pub fn normalized_entropy(&self) -> Option { if !self.ready() { return None;