OBIKmers/obikmer
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 27 Wiki Activity

docs: add obikmer user guide and MkDocs build configuration

Browse Source 
Introduces a comprehensive documentation set covering theoretical foundations, CLI usage, installation, and system architecture. Adds MkDocs configuration and Makefile targets to generate, serve with live reload, and clean the documentation site. Includes citation styles and bibliography files for academic references.
This commit is contained in:
Eric Coissac
2026-08-13 17:19:01 +02:00
parent 0c86ea0385
commit 6acafa7f2c
28 changed files with 1497 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
UserDocMD/usage/distance.md
+97
View File
@@ -0,0 +1,97 @@
# distance
Compute pairwise distance matrices between the genomes stored in an index, optionally build trees (NJ/UPGMA), and optionally derive a central-position SNP model with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE).
```bash
obikmer distance INDEX [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory |
## Distance matrix
| Option | Default | Description |
|---|---|---|
| `--metric` | `jaccard` | One of `jaccard`, `mash`, `hamming`, `bray-curtis`, `relfreq-bray-curtis`, `euclidean`, `relfreq-euclidean`, `hellinger`, `hellinger-euclidean` |
| `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for Jaccard/Mash on a count index |
| `--shared-kmers` | off | Also write the shared-kmer count matrix |
| `--nj` | off | Compute and write a Neighbor-Joining tree (Newick) |
| `--upgma` | off | Compute and write a UPGMA tree (Newick) |
| `-o, --output` | none (stdout) | Output file prefix; without it, the distance matrix is printed to stdout as CSV |
`hamming` requires a presence/absence index. All other metrics work on either index type; on a presence index, `jaccard`/`mash`/`hamming` are the only ones available.
### Metric definitions
- **jaccard**: $D = 1 - \dfrac{|A \cap B|}{|A \cup B|}$ over the sets of kmers present in each genome.
- **mash**: derived from the Jaccard distance via $D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)$ where $J = 1 - D_{\text{jaccard}}$ and $k$ is the index's kmer size; clamped to 1.0 when $J \le 0$.
- **hamming**: number of kmer positions where presence differs between the two genomes (presence index only, not normalized): $D = \sum_i \mathbb{1}[a_i \ne b_i]$.
- **bray-curtis**: $D = 1 - \dfrac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}$ on raw per-kmer counts.
- **relfreq-bray-curtis**: the same formula computed on per-genome relative frequencies $p_i = c_i / \sum_j c_j$ instead of raw counts.
- **euclidean**: $D = \sqrt{\sum_i (c_i^A - c_i^B)^2}$ on raw counts.
- **relfreq-euclidean**: the same formula on relative frequencies.
- **hellinger**: $D = \dfrac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}$ on relative frequencies, bounded in $[0, 1]$.
- **hellinger-euclidean**: the unnormalized variant, $D = \sqrt{2} \times D_{\text{hellinger}}$.
## Central-position SNP model
This is a separate operation from the distance-matrix computation above: if any option below is used, no `--metric` matrix is computed in the same invocation.
A **family** is the set of up to 4 kmers that share identical flanking sequence and differ only at the exact central base. Because $k$ is odd, the central position is well defined and maps to itself under reverse complementation. All computations below first require building the **sibling annex**, an index-wide record of which of the 4 possible central bases are observed at each family, across every genome.
| Option | Default | Description |
|---|---|---|
| `--sibling-annex` | off | Build the sibling presence-mask annex (prerequisite for every option below) |
| `--exclude-genome LABEL` | none | Exclude a genome (repeatable) from every SNP/Sankoff/export computation below |
| `--sibling-stats` | off | Write the family-size (sibling count) distribution, per genome and globally |
| `--raw-snp-distance` | off | Write the single-copy central-SNP p-distance matrix |
| `--raw-snp-counts` | off | Write per-pair diagnostic counts (n_snp, n_shared, n_eligible) instead of a matrix |
| `--snp` | off | Write a SNP-only pseudo-alignment in FASTA, IUPAC-coded |
### Locus eligibility
A family is eligible for a genome pair $(i, j)$ only if genome $i$ carries exactly one of the family's observed forms (single-copy, unambiguous) and genome $j$ also carries exactly one. A genome carrying more than one form at a locus makes that locus ineligible for any pair involving it.
`--raw-snp-distance` tallies, over every eligible locus of every genome pair, $n_{\text{snp}}$ (the two genomes' single forms differ) versus $n_{\text{shared}}$ (they agree — this includes invariant families). The output ratio is $\hat{p} = \dfrac{n_{\text{snp}}}{n_{\text{snp}} + n_{\text{shared}}}$.
`--snp` restricts itself to *variable* families (family size $\ge 2$) and writes one FASTA record per genome, one column per family, IUPAC-coded from each genome's presence mask at that family (a single form → the plain base; several forms → the matching IUPAC ambiguity code; no form → `-`).
`--exclude-genome` removes a genome from these computations, re-checking column variability among the remaining genomes so that a column made monomorphic by the exclusion is dropped rather than kept artificially. It does not affect the `--metric` distance-matrix path.
## Sankoff calibration and phylogenetic exports
| Option | Default | Description |
|---|---|---|
| `--sankoff` | off | Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment |
| `--sankoff-ratio-ceiling` | `0.5` | Exclude genome pairs whose raw SNP ratio exceeds this value from the calibration |
| `--tnt` | off | Also write a TNT script (implies `--sankoff`) |
| `--phyg` | off | Also write PhyG input files (implies `--sankoff`) |
| `--iqtree` | off | Also write an IQ-TREE custom model and alignment (implies `--sankoff`) |
| `--sankoff-cost-scale` | `100` | Integer scaling factor applied to costs before rounding (required by TNT/PhyG's integer-only cost commands) |
### The 16-state model
Each family is treated as a character with 16 possible states: one per subset of the 4 possible central bases actually observed (including the empty subset). Calibration combines two tallies, both restricted to genome pairs at or below `--sankoff-ratio-ceiling`:
- a $5 \times 5$ transition matrix over family cardinality (0–4 observed forms) between paired genomes, and
- a $4 \times 4$ base-substitution transition matrix from unambiguous single-copy loci,
which are combined into a row-normalized $16 \times 16$ transition probability matrix $P$, converted to a symmetric cost matrix via $\text{cost}(a,b) = -\ln P(a,b)$.
`--sankoff` alone writes the cost matrix, the calibration parameters, and a pseudo-alignment recoded so the empty state uses the symbol `0` (never a gap character, to avoid ambiguity with external tools' own gap semantics). It does not run any external tool.
### Exports
All three exports reuse the `--sankoff` calibrated matrix and pseudo-alignment, recoded for the target tool:
- **`--tnt`**: a self-contained TNT script (alignment recoded to TNT's fixed 16-symbol alphabet, integer-scaled cost matrix re-closed to a metric, a default search block).
- **`--phyg`**: a custom cost-matrix file plus a PhyG script reusing the `--sankoff` alignment directly.
- **`--iqtree`**: a custom substitution-model file (exchangeability matrix recovered as $R(a,b) = e^{-\text{cost}(a,b)}$, plus empirical state frequencies) and a matching alignment, for maximum-likelihood inference with real branch lengths (unlike the parsimony step-counts from TNT/PhyG). Only states actually occurring in the alignment are kept and compactly renumbered.
## Output files
With `-o/--output PREFIX`, the relevant subset of the following files is written: `<prefix>_dist.csv`, `<prefix>_shared.csv`, `<prefix>_nj.nwk`, `<prefix>_upgma.nwk`, `<prefix>_siblings.csv`, `<prefix>_rawsnp.csv`, `<prefix>_rawsnp_counts.csv`, `<prefix>_snp.fasta`, `<prefix>_sankoff_matrix.csv`, `<prefix>_sankoff_params.yaml`, `<prefix>_sankoff.fasta`, `<prefix>_sankoff.tnt`, `<prefix>_sankoff.tcm`, `<prefix>_sankoff.pg`, `<prefix>_iqtree.model`, `<prefix>_iqtree.fasta`. Without `-o`, only the plain `--metric` distance matrix is produced, on stdout.