The 16-state Sankoff model

Each family (see Central-position SNP model) is treated as a character with 16 possible states: one per subset of the 4 possible central bases, including the empty subset (no observed form). A genome pair's calibration reduces to two separately-estimated first-order Markov models, both restricted to genome pairs below a configurable SNP-ratio ceiling (excluding saturated pairs from skewing the calibration):

Cardinality transitions (5 \times 5, row-stochastic): the probability that a family observed with i forms in one genome (i \in \{0,\dots,4\}) is observed with j forms in the other.
Composition transitions (4 \times 4, row-stochastic): the probability that an unambiguous, single-copy base observed as one of A/C/G/T in one genome is observed as another base in the other genome.

Combining into a 16×16 cost matrix

The two matrices are not combined as a plain tensor product. For a pair of states (A, B) — each a subset of \{A,C,G,T\} — let shared = A \cap B, lost = A \setminus B, gained = B \setminus A. The log-probability of the transition A \to B is:

\log P(A,B) = \underbrace{\log P_{\text{card}}(|A|,|B|)}_{\text{cardinality term, omitted when loss events are treated as free}} \;+\; \underbrace{\sum_{i \,\in\, \text{shared}} \log P_{\text{comp}}(i,i)}_{\text{bases conserved on both sides}} \;-\; \underbrace{\text{best\_pairing\_cost}(\text{lost}, \text{gained})}_{\text{substitutions among the differing bases}}

best_pairing_cost enumerates every injection pairing elements of lost with elements of gained and keeps the one minimizing -\sum \log P_{\text{comp}} over the pairs — the discrete analogue of "prefer a substitution over an independent loss and gain": if a base is lost from one side and a different base is gained on the other, that is scored as a single substitution between them (cheaper, under any sane calibration, than treating them as two unrelated cardinality-changing events) whenever such a pairing is possible.

The resulting 16 \times 16 log-probability matrix is row-normalized in log-space (log-sum-exp, not a direct exp()/divide, for numerical stability), giving a genuine row-stochastic transition matrix P. The cost matrix is then \text{cost}(A,B) = -\ln P(A,B), symmetrized by simple averaging:

\text{sym}(A,B) = \frac{\text{cost}(A,B) + \text{cost}(B,A)}{2}

This calibrated cost matrix is what every downstream export (for TNT, PhyG, IQ-TREE) reuses, recoded into each tool's own format.

See phylo for how to run the calibration and produce these exports.

Wiki sidebar

Theory

Kmer indexing

Phylogeny

Kmer-based

Whole k-mer set distances