4123 lines
186 KiB
HTML
4123 lines
186 KiB
HTML
|
||
<!doctype html>
|
||
<html lang="en" class="no-js">
|
||
<head>
|
||
|
||
<meta charset="utf-8">
|
||
<meta name="viewport" content="width=device-width,initial-scale=1">
|
||
|
||
|
||
|
||
|
||
<link rel="prev" href="../indexing/">
|
||
|
||
|
||
<link rel="next" href="../../implementation/superkmer/">
|
||
|
||
|
||
|
||
|
||
|
||
<link rel="icon" href="../../assets/images/favicon.png">
|
||
<meta name="generator" content="mkdocs-1.6.1, mkdocs-material-9.7.6">
|
||
|
||
|
||
|
||
<title>Central-position SNP distance (discussion) - obikmer</title>
|
||
|
||
|
||
|
||
<link rel="stylesheet" href="../../assets/stylesheets/main.484c7ddc.min.css">
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<link rel="preconnect" href="https://fonts.gstatic.com" crossorigin>
|
||
<link rel="stylesheet" href="https://fonts.googleapis.com/css?family=Roboto:300,300i,400,400i,700,700i%7CRoboto+Mono:400,400i,700,700i&display=fallback">
|
||
<style>:root{--md-text-font:"Roboto";--md-code-font:"Roboto Mono"}</style>
|
||
|
||
|
||
|
||
<script>__md_scope=new URL("../..",location),__md_hash=e=>[...e].reduce(((e,_)=>(e<<5)-e+_.charCodeAt(0)),0),__md_get=(e,_=localStorage,t=__md_scope)=>JSON.parse(_.getItem(t.pathname+"."+e)),__md_set=(e,_,t=localStorage,a=__md_scope)=>{try{t.setItem(a.pathname+"."+e,JSON.stringify(_))}catch(e){}}</script>
|
||
|
||
|
||
|
||
|
||
|
||
</head>
|
||
|
||
|
||
<body dir="ltr">
|
||
|
||
|
||
<input class="md-toggle" data-md-toggle="drawer" type="checkbox" id="__drawer" autocomplete="off">
|
||
<input class="md-toggle" data-md-toggle="search" type="checkbox" id="__search" autocomplete="off">
|
||
<label class="md-overlay" for="__drawer"></label>
|
||
<div data-md-component="skip">
|
||
|
||
|
||
<a href="#central-position-snp-distance-discussion" class="md-skip">
|
||
Skip to content
|
||
</a>
|
||
|
||
</div>
|
||
<div data-md-component="announce">
|
||
|
||
</div>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<header class="md-header md-header--shadow" data-md-component="header">
|
||
<nav class="md-header__inner md-grid" aria-label="Header">
|
||
<a href="../.." title="obikmer" class="md-header__button md-logo" aria-label="obikmer" data-md-component="logo">
|
||
|
||
|
||
<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 24 24"><path d="M12 8a3 3 0 0 0 3-3 3 3 0 0 0-3-3 3 3 0 0 0-3 3 3 3 0 0 0 3 3m0 3.54C9.64 9.35 6.5 8 3 8v11c3.5 0 6.64 1.35 9 3.54 2.36-2.19 5.5-3.54 9-3.54V8c-3.5 0-6.64 1.35-9 3.54"/></svg>
|
||
|
||
</a>
|
||
<label class="md-header__button md-icon" for="__drawer">
|
||
|
||
<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 24 24"><path d="M3 6h18v2H3zm0 5h18v2H3zm0 5h18v2H3z"/></svg>
|
||
</label>
|
||
<div class="md-header__title" data-md-component="header-title">
|
||
<div class="md-header__ellipsis">
|
||
<div class="md-header__topic">
|
||
<span class="md-ellipsis">
|
||
obikmer
|
||
</span>
|
||
</div>
|
||
<div class="md-header__topic" data-md-component="header-topic">
|
||
<span class="md-ellipsis">
|
||
|
||
Central-position SNP distance (discussion)
|
||
|
||
</span>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
|
||
<script>var palette=__md_get("__palette");if(palette&&palette.color){if("(prefers-color-scheme)"===palette.color.media){var media=matchMedia("(prefers-color-scheme: light)"),input=document.querySelector(media.matches?"[data-md-color-media='(prefers-color-scheme: light)']":"[data-md-color-media='(prefers-color-scheme: dark)']");palette.color.media=input.getAttribute("data-md-color-media"),palette.color.scheme=input.getAttribute("data-md-color-scheme"),palette.color.primary=input.getAttribute("data-md-color-primary"),palette.color.accent=input.getAttribute("data-md-color-accent")}for(var[key,value]of Object.entries(palette.color))document.body.setAttribute("data-md-color-"+key,value)}</script>
|
||
|
||
|
||
|
||
|
||
</nav>
|
||
|
||
</header>
|
||
|
||
<div class="md-container" data-md-component="container">
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<main class="md-main" data-md-component="main">
|
||
<div class="md-main__inner md-grid">
|
||
|
||
|
||
|
||
<div class="md-sidebar md-sidebar--primary" data-md-component="sidebar" data-md-type="navigation" >
|
||
<div class="md-sidebar__scrollwrap">
|
||
<div class="md-sidebar__inner">
|
||
|
||
|
||
|
||
|
||
<nav class="md-nav md-nav--primary" aria-label="Navigation" data-md-level="0">
|
||
<label class="md-nav__title" for="__drawer">
|
||
<a href="../.." title="obikmer" class="md-nav__button md-logo" aria-label="obikmer" data-md-component="logo">
|
||
|
||
|
||
<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 24 24"><path d="M12 8a3 3 0 0 0 3-3 3 3 0 0 0-3-3 3 3 0 0 0-3 3 3 3 0 0 0 3 3m0 3.54C9.64 9.35 6.5 8 3 8v11c3.5 0 6.64 1.35 9 3.54 2.36-2.19 5.5-3.54 9-3.54V8c-3.5 0-6.64 1.35-9 3.54"/></svg>
|
||
|
||
</a>
|
||
obikmer
|
||
</label>
|
||
|
||
<ul class="md-nav__list" data-md-scrollfix>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../.." class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Home
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../installation/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Installation
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item md-nav__item--active md-nav__item--nested">
|
||
|
||
|
||
|
||
<input class="md-nav__toggle md-toggle " type="checkbox" id="__nav_3" checked>
|
||
|
||
|
||
<label class="md-nav__link" for="__nav_3" id="__nav_3_label" tabindex="0">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Theory
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
<span class="md-nav__icon md-icon"></span>
|
||
</label>
|
||
|
||
<nav class="md-nav" data-md-level="1" aria-labelledby="__nav_3_label" aria-expanded="true">
|
||
<label class="md-nav__title" for="__nav_3">
|
||
<span class="md-nav__icon md-icon"></span>
|
||
|
||
|
||
Theory
|
||
|
||
|
||
</label>
|
||
<ul class="md-nav__list" data-md-scrollfix>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../kmers/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Kmers and super-kmers
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../encoding/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
DNA encoding
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../entropy/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Entropy filter
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../minimizer/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Minimizer selection
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../indexing/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Partitioning architecture
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item md-nav__item--active">
|
||
|
||
<input class="md-nav__toggle md-toggle" type="checkbox" id="__toc">
|
||
|
||
|
||
|
||
|
||
|
||
<label class="md-nav__link md-nav__link--active" for="__toc">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Central-position SNP distance (discussion)
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
<span class="md-nav__icon md-icon"></span>
|
||
</label>
|
||
|
||
<a href="./" class="md-nav__link md-nav__link--active">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Central-position SNP distance (discussion)
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
|
||
|
||
|
||
<nav class="md-nav md-nav--secondary" aria-label="Table of contents">
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<label class="md-nav__title" for="__toc">
|
||
<span class="md-nav__icon md-icon"></span>
|
||
Table of contents
|
||
</label>
|
||
<ul class="md-nav__list" data-md-component="toc" data-md-scrollfix>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#motivation" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Motivation
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#statistic-and-correspondence-with-shared" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Statistic and correspondence with shared
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Statistic and correspondence with shared">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#definitions-family-and-the-canonical-form-of-a-family" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Definitions: family, and the canonical form of a family
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#locus-eligibility-raw-definition-vs-paralogy-filter" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Locus eligibility: raw definition vs. paralogy filter
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#multi-genome-framing-family-as-pseudo-alignment-column" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Multi-genome framing: family as pseudo-alignment column
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#context-detectability-and-a-3-way-ordinal-distance-per-pair" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Context, detectability, and a 3-way ordinal distance per pair
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#sankoff-parsimony-as-the-resolution-of-the-16-state-model-problem" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Sankoff parsimony as the resolution of the 16-state model problem
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Sankoff parsimony as the resolution of the 16-state model problem">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#a-concrete-sankoff-cost-matrix-for-the-16-state-alphabet" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
A concrete Sankoff cost matrix for the 16-state alphabet
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#experiment-tnt-run-on-real-data-2026-08-11" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Experiment: TNT run on real data (2026-08-11)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#native-sankoff-tnt-phyg-export-2026-08-12-superseding-the-external-scripts-above" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Native --sankoff --tnt/--phyg export (2026-08-12), superseding the external scripts above
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#next-direction-genuine-ml-branch-lengths-not-parsimony-open-2026-08-12" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Next direction: genuine ML branch lengths, not parsimony (open, 2026-08-12)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#r-implemented-exclude-genome-added-rogue-taxon-test-negative-2026-08-12" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
R/π implemented; --exclude-genome added; rogue-taxon test negative (2026-08-12)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#r-via-exp-cost-is-wrong-for-a-ctmc-cardinalitycomposition-decomposition-open-2026-08-12" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
R via exp(-cost) is wrong for a CTMC; cardinality/composition decomposition (open, 2026-08-12)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#implemented-2026-08-12-pairwise_cost_matrix-replaces-build_cost_matrix-entirely" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Implemented (2026-08-12): pairwise_cost_matrix replaces build_cost_matrix entirely
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#two-consistency-bugs-found-and-fixed-post-implementation-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Two consistency bugs found and fixed post-implementation (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#heterozygosity-ploidy-and-consensus-assembly-inputs" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Heterozygosity, ploidy, and consensus-assembly inputs
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#sufficient-statistic-4x4-base-pair-tally" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Sufficient statistic: 4x4 base-pair tally
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#biases-properties-of-the-estimator-not-defects" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Biases (properties of the estimator, not defects)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#implementation-avoid-materializing-a-de-bruijn-graph" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Implementation: avoid materializing a de Bruijn graph
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#implementation-sequential-per-partition-sweep-no-scratch-no-graph" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Implementation: sequential per-partition sweep (no scratch, no graph)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#cheaper-subsampling" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Cheaper: subsampling
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#recommendation" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Recommendation
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#detailed-implementation-plan" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Detailed implementation plan
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Detailed implementation plan">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-0-new-low-level-primitives-obikseq" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 0 — new low-level primitives (obikseq)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-1-the-tally-accumulator-obikindex" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 1 — the tally accumulator (obikindex)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-2-the-sweep-obikindex-new-snprs" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 2 — the sweep (obikindex, new snp.rs)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#open-problem-unresolved-session-end-not-yet-fully-convinced" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Open problem (unresolved, session end — not yet fully convinced)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-2b-sibling-count-minorant-annex-consolidated-plan" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 2b — sibling-count / minorant annex (consolidated plan)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-3-finalisation-obikindex" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 3 — finalisation (obikindex)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-4-surfacing-obikindex-obikmer-cli" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 4 — surfacing (obikindex + obikmer CLI)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-5-subsampling-flag" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 5 — subsampling flag
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#testing" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Testing
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#suggested-phasing" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Suggested phasing
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#locus-dropout-under-incomplete-coverage-absence-as-missing-data-not-a-penalized-state-2026-08-13-discussion" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Locus dropout under incomplete coverage: absence as missing data, not a penalized state (2026-08-13, discussion)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Locus dropout under incomplete coverage: absence as missing data, not a penalized state (2026-08-13, discussion)">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#-free-loss-implemented-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
--free-loss implemented (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#tnt-confirmed-2026-08-13-iq-tree-asc-constant-column-fix-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
TNT confirmed (2026-08-13); IQ-TREE +ASC constant-column fix (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#-family-overlap-min-shared-family-and-a-second-free_loss-bug-cardinality-transition-cost-not-zeroed-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
--family-overlap/--min-shared-family, and a second free_loss bug: cardinality-transition cost not zeroed (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#references" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
References
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
|
||
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item md-nav__item--nested">
|
||
|
||
|
||
|
||
<input class="md-nav__toggle md-toggle " type="checkbox" id="__nav_4" >
|
||
|
||
|
||
<label class="md-nav__link" for="__nav_4" id="__nav_4_label" tabindex="0">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Implementation
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
<span class="md-nav__icon md-icon"></span>
|
||
</label>
|
||
|
||
<nav class="md-nav" data-md-level="1" aria-labelledby="__nav_4_label" aria-expanded="false">
|
||
<label class="md-nav__title" for="__nav_4">
|
||
<span class="md-nav__icon md-icon"></span>
|
||
|
||
|
||
Implementation
|
||
|
||
|
||
</label>
|
||
<ul class="md-nav__list" data-md-scrollfix>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/superkmer/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
SuperKmer
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/kmer/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Kmer
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/chunkreader/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Chunk reader
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/pipeline/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Construction pipeline
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/obipipeline/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
obipipeline library
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/storage/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
On-disk storage
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/mphf/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
MPHF selection
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/unitig_evidence/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Unitig evidence encoding
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/evidence_elimination/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Evidence elimination (discussion)
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/obilayeredmap/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
obilayeredmap crate
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/persistent_compact_int_vec/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
PersistentCompactIntVec
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/persistent_bit_vec/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
PersistentBitVec
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/merge/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Merge command
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/merge_parallelism/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Merge parallelism & memory
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/filtering/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Kmer filtering
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/select/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Select command
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../implementation/obitaxonomy/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
obitaxonomy crate
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item md-nav__item--nested">
|
||
|
||
|
||
|
||
<input class="md-nav__toggle md-toggle " type="checkbox" id="__nav_5" >
|
||
|
||
|
||
<label class="md-nav__link" for="__nav_5" id="__nav_5_label" tabindex="0">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Architecture
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
<span class="md-nav__icon md-icon"></span>
|
||
</label>
|
||
|
||
<nav class="md-nav" data-md-level="1" aria-labelledby="__nav_5_label" aria-expanded="false">
|
||
<label class="md-nav__title" for="__nav_5">
|
||
<span class="md-nav__icon md-icon"></span>
|
||
|
||
|
||
Architecture
|
||
|
||
|
||
</label>
|
||
<ul class="md-nav__list" data-md-scrollfix>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../architecture/sequences/invariant/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Sequences
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../architecture/index_architecture/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Kmer index
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../architecture/siblings/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
Sibling annex (discussion)
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../architecture/numa_worker_pools/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
NUMA-aware worker pools
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<li class="md-nav__item">
|
||
<a href="../../architecture/numa_partition_runner/" class="md-nav__link">
|
||
|
||
|
||
|
||
<span class="md-ellipsis">
|
||
|
||
|
||
NUMA-aware partition runner
|
||
|
||
|
||
|
||
</span>
|
||
|
||
|
||
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
|
||
|
||
</ul>
|
||
</nav>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
|
||
|
||
<div class="md-sidebar md-sidebar--secondary" data-md-component="sidebar" data-md-type="toc" >
|
||
<div class="md-sidebar__scrollwrap">
|
||
<div class="md-sidebar__inner">
|
||
|
||
|
||
<nav class="md-nav md-nav--secondary" aria-label="Table of contents">
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<label class="md-nav__title" for="__toc">
|
||
<span class="md-nav__icon md-icon"></span>
|
||
Table of contents
|
||
</label>
|
||
<ul class="md-nav__list" data-md-component="toc" data-md-scrollfix>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#motivation" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Motivation
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#statistic-and-correspondence-with-shared" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Statistic and correspondence with shared
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Statistic and correspondence with shared">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#definitions-family-and-the-canonical-form-of-a-family" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Definitions: family, and the canonical form of a family
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#locus-eligibility-raw-definition-vs-paralogy-filter" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Locus eligibility: raw definition vs. paralogy filter
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#multi-genome-framing-family-as-pseudo-alignment-column" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Multi-genome framing: family as pseudo-alignment column
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#context-detectability-and-a-3-way-ordinal-distance-per-pair" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Context, detectability, and a 3-way ordinal distance per pair
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#sankoff-parsimony-as-the-resolution-of-the-16-state-model-problem" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Sankoff parsimony as the resolution of the 16-state model problem
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Sankoff parsimony as the resolution of the 16-state model problem">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#a-concrete-sankoff-cost-matrix-for-the-16-state-alphabet" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
A concrete Sankoff cost matrix for the 16-state alphabet
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#experiment-tnt-run-on-real-data-2026-08-11" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Experiment: TNT run on real data (2026-08-11)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#native-sankoff-tnt-phyg-export-2026-08-12-superseding-the-external-scripts-above" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Native --sankoff --tnt/--phyg export (2026-08-12), superseding the external scripts above
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#next-direction-genuine-ml-branch-lengths-not-parsimony-open-2026-08-12" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Next direction: genuine ML branch lengths, not parsimony (open, 2026-08-12)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#r-implemented-exclude-genome-added-rogue-taxon-test-negative-2026-08-12" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
R/π implemented; --exclude-genome added; rogue-taxon test negative (2026-08-12)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#r-via-exp-cost-is-wrong-for-a-ctmc-cardinalitycomposition-decomposition-open-2026-08-12" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
R via exp(-cost) is wrong for a CTMC; cardinality/composition decomposition (open, 2026-08-12)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#implemented-2026-08-12-pairwise_cost_matrix-replaces-build_cost_matrix-entirely" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Implemented (2026-08-12): pairwise_cost_matrix replaces build_cost_matrix entirely
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#two-consistency-bugs-found-and-fixed-post-implementation-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Two consistency bugs found and fixed post-implementation (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#heterozygosity-ploidy-and-consensus-assembly-inputs" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Heterozygosity, ploidy, and consensus-assembly inputs
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#sufficient-statistic-4x4-base-pair-tally" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Sufficient statistic: 4x4 base-pair tally
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#biases-properties-of-the-estimator-not-defects" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Biases (properties of the estimator, not defects)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#implementation-avoid-materializing-a-de-bruijn-graph" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Implementation: avoid materializing a de Bruijn graph
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#implementation-sequential-per-partition-sweep-no-scratch-no-graph" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Implementation: sequential per-partition sweep (no scratch, no graph)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#cheaper-subsampling" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Cheaper: subsampling
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#recommendation" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Recommendation
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#detailed-implementation-plan" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Detailed implementation plan
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Detailed implementation plan">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-0-new-low-level-primitives-obikseq" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 0 — new low-level primitives (obikseq)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-1-the-tally-accumulator-obikindex" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 1 — the tally accumulator (obikindex)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-2-the-sweep-obikindex-new-snprs" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 2 — the sweep (obikindex, new snp.rs)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#open-problem-unresolved-session-end-not-yet-fully-convinced" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Open problem (unresolved, session end — not yet fully convinced)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-2b-sibling-count-minorant-annex-consolidated-plan" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 2b — sibling-count / minorant annex (consolidated plan)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-3-finalisation-obikindex" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 3 — finalisation (obikindex)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-4-surfacing-obikindex-obikmer-cli" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 4 — surfacing (obikindex + obikmer CLI)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#step-5-subsampling-flag" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Step 5 — subsampling flag
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#testing" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Testing
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#suggested-phasing" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Suggested phasing
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#locus-dropout-under-incomplete-coverage-absence-as-missing-data-not-a-penalized-state-2026-08-13-discussion" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
Locus dropout under incomplete coverage: absence as missing data, not a penalized state (2026-08-13, discussion)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
<nav class="md-nav" aria-label="Locus dropout under incomplete coverage: absence as missing data, not a penalized state (2026-08-13, discussion)">
|
||
<ul class="md-nav__list">
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#-free-loss-implemented-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
--free-loss implemented (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#tnt-confirmed-2026-08-13-iq-tree-asc-constant-column-fix-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
TNT confirmed (2026-08-13); IQ-TREE +ASC constant-column fix (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#-family-overlap-min-shared-family-and-a-second-free_loss-bug-cardinality-transition-cost-not-zeroed-2026-08-13" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
--family-overlap/--min-shared-family, and a second free_loss bug: cardinality-transition cost not zeroed (2026-08-13)
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
</nav>
|
||
|
||
</li>
|
||
|
||
<li class="md-nav__item">
|
||
<a href="#references" class="md-nav__link">
|
||
<span class="md-ellipsis">
|
||
|
||
References
|
||
|
||
</span>
|
||
</a>
|
||
|
||
</li>
|
||
|
||
</ul>
|
||
|
||
</nav>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
|
||
|
||
<div class="md-content" data-md-component="content">
|
||
|
||
<article class="md-content__inner md-typeset">
|
||
|
||
|
||
|
||
|
||
|
||
<h1 id="central-position-snp-distance-discussion">Central-position SNP distance (discussion)</h1>
|
||
<p>Not implemented. Design discussion for a substitution-rate estimator that
|
||
observes SNPs directly from paired-genome k-mer comparison, as an alternative
|
||
to Mash's Poisson-Jaccard inference (see <a href="../../implementation/obicompactvec/">obicompactvec</a>
|
||
for the implemented Jaccard/Mash distances).</p>
|
||
<h2 id="motivation">Motivation</h2>
|
||
<p><strong>Primary intent: restrict the comparison to what is actually comparable.</strong>
|
||
Mash's Jaccard is computed over the <strong>union</strong> of both genomes' k-mer content:
|
||
anything not identically shared is folded into a single undifferentiated
|
||
mass, whether the cause is a point substitution, a genuinely absent
|
||
homologous region (lineage-specific content, gene-family expansion, HGT,
|
||
genome-size asymmetry), or a diverged paralogous copy. The model then
|
||
back-infers a single mutation rate from that mass, silently attributing
|
||
non-homology to mutation. The central-SNP approach instead conditions every
|
||
comparison on local, positive evidence of homology: a locus only enters the
|
||
statistic if its <code>2m</code> flanking bases (<code>m = (k-1)/2</code>) are found intact in
|
||
<em>both</em> genomes — genuinely absent or non-homologous content is excluded from
|
||
the comparison entirely (neither numerator nor denominator), rather than
|
||
silently counted as divergence. This is a conditioning on comparability, not
|
||
just a richer summary statistic — see "Statistic and correspondence with
|
||
<code>shared</code>" below for how it plays out against genome-size asymmetry and
|
||
diverged gene families, and "Heterozygosity, ploidy, and consensus-assembly
|
||
inputs" for the corresponding paralogy/heterozygosity filter.</p>
|
||
<p><strong>Secondary benefit: access to the substitution's nature.</strong> Because the
|
||
central base of an odd-k window is directly observable once the flanks are
|
||
confirmed conserved, this also yields more than a rate — the
|
||
transition/transversion split — enabling classical corrected distances
|
||
(Jukes-Cantor, Kimura 2-parameter, LogDet) that a single Jaccard scalar
|
||
cannot support.</p>
|
||
<h2 id="statistic-and-correspondence-with-shared">Statistic and correspondence with <code>shared</code></h2>
|
||
<p>A genomic position <code>p</code> is covered by <code>k</code> overlapping k-mer windows. Requiring
|
||
the substitution to sit at the window's <strong>center</strong> makes exactly one window
|
||
per SNP eligible — a 1:1 correspondence between SNP and center-neighbor k-mer
|
||
pair, avoiding the ~k-fold overcount of an any-position neighbor search.</p>
|
||
<p>A locus with a fully conserved <code>k</code>-window (flanks <strong>and</strong> center) is an
|
||
exact-shared k-mer at that locus; a locus with conserved flanks but a
|
||
substituted center is a "central SNP". Both count each locus exactly once, in
|
||
matching units:</p>
|
||
<div class="highlight"><pre><span></span><code>p_hat[i,j] = SNP[i,j] / (SNP[i,j] + shared[i,j])
|
||
</code></pre></div>
|
||
<p><code>p_hat</code> is <code>P(center substituted | 2m flanks conserved)</code>. <code>shared[i,j]</code> here
|
||
is <strong>not</strong> the general-purpose <code>shared_kmers</code> matrix used by Jaccard/Mash
|
||
(<code>--shared-kmers</code>, <code>BitPartials::partial_jaccard</code> /
|
||
<code>CountPartials::partial_threshold_jaccard</code>) — that matrix counts raw k-mer
|
||
identity with no per-genome copy-number constraint, whereas <code>p_hat</code>'s
|
||
denominator applies the eligibility rule defined below (raw or
|
||
paralogy-filtered). Both <code>SNP</code> and <code>shared</code> are accumulated by the same
|
||
sweep, from the same per-locus candidate set (source k-mer + 3 variants),
|
||
under the same eligibility rule — see "Locus eligibility" below, and
|
||
"Heterozygosity, ploidy, and consensus-assembly inputs" for why the
|
||
copy-number constraint matters and what it costs.</p>
|
||
<p><strong>Canonical invariance</strong>: for odd k, the central position maps to itself under
|
||
reverse-complement (<code>m -> k-1-m = m</code>, base complemented). A transition maps to
|
||
a transition, a transversion to a transversion — the transition/transversion
|
||
split is well-defined in canonical space.</p>
|
||
<h3 id="definitions-family-and-the-canonical-form-of-a-family">Definitions: family, and the canonical form of a family</h3>
|
||
<p><strong>Family.</strong> The family of a k-mer <code>x</code> is the set of (up to) 4 k-mers sharing
|
||
<code>x</code>'s <code>2m</code> flanking bases, differing only at the central base <code>m</code>. Membership
|
||
is a property of the flank pattern, not of <code>x</code> itself: any of the 4 possible
|
||
central substitutions belongs to the same family.</p>
|
||
<p><strong><code>central_canonical_neighbors()</code></strong> (<code>obikseq</code>, <code>CanonicalKmerOf::central_canonical_neighbors</code>)
|
||
generates all 4 members from any one of them (observed or not), each
|
||
independently canonicalised (<code>.canonical()</code>, i.e. <code>min(kmer, revcomp(kmer))</code>).
|
||
This independent canonicalisation is necessary because a central substitution
|
||
can flip which orientation is lexicographically smaller — two members of the
|
||
same family can end up canonicalised in <em>different</em> orientations. Despite
|
||
that, the <strong>set</strong> of 4 resulting canonical k-mers is invariant: calling
|
||
<code>central_canonical_neighbors()</code> on any member of a family — present in the
|
||
index or not — yields the same 4 values. This is relied upon throughout the
|
||
rest of this document.</p>
|
||
<p><strong>Canonical form of a family.</strong> Because orientation can differ member to
|
||
member, "which of the 4 is the reference" cannot be defined relative to
|
||
<em>whichever member happened to be visited first</em>, nor relative to the
|
||
minorant (see below) — both are data-dependent (they depend on what is
|
||
actually observed), so using either as the reference would make the
|
||
reference itself vary depending on what happens to be present in a given
|
||
index. Instead: <strong>the canonical form of a family is, by definition, the
|
||
member whose own central base — read in its own already-canonical
|
||
orientation — is <code>A</code>.</strong> This is well-defined for every family, computed
|
||
purely from the flank pattern, whether or not that specific member (or any
|
||
member at all) is actually observed anywhere in the index. Concretely: call
|
||
<code>central_canonical_neighbors()</code> on any member (observed or not) to get the
|
||
family's 4 canonical forms; the one among them whose own centre nucleotide is
|
||
<code>A</code> is the family's canonical form. The other 3 (<code>C</code>, <code>G</code>, <code>T</code>) are labelled
|
||
relative to <em>that</em> fixed reference, not relative to the calling member's own
|
||
orientation.</p>
|
||
<p><strong>Consequence for the minorant.</strong> With this fixed A-referenced labelling,
|
||
<code>minorant</code> (the smallest raw encoding among the family's <em>observed</em> members,
|
||
introduced further below) becomes directly computable rather than needing to
|
||
be tracked as extra state: regenerate the family's 4 canonical forms from
|
||
any member's own k-mer (cheap, no lookup), compare the raw encodings of
|
||
whichever are marked present, and take the smallest. No separate stored bit
|
||
is required — see Step 2b below, where this replaces the earlier
|
||
minorant-bit design.</p>
|
||
<h2 id="locus-eligibility-raw-definition-vs-paralogy-filter">Locus eligibility: raw definition vs. paralogy filter</h2>
|
||
<p>For each k-mer <code>x</code> observed in genome A (source, one MPHF slot; the 3
|
||
central-position variants generated as in the sweep below): check whether
|
||
A's locus (flanks fixed) is resolvable in genome B under one of the 4
|
||
central forms.</p>
|
||
<p><strong>Raw / no model.</strong> The locus counts in the denominator iff at least one of
|
||
the 4 forms is present in B; it counts in the numerator iff the form found in
|
||
B differs from A's own. No constraint on A's or B's own copy number at this
|
||
locus. Open question, not resolved: what if <strong>more than one</strong> of the 4 forms
|
||
is present in B simultaneously (ambiguous target — count once arbitrarily,
|
||
count all, or drop)? The stringent filter below sidesteps the question by
|
||
construction rather than answering it.</p>
|
||
<p><strong>Stringent / paralogy-aware.</strong> The locus counts only if exactly one of the
|
||
4 forms is present in A <strong>and</strong> exactly one is present in B (<code>count == 1</code> at
|
||
that slot too, when a count index is available, to also exclude same-allele
|
||
duplicates that presence alone cannot see). This drops the raw definition's
|
||
ambiguous-B case automatically, at the cost of also dropping heterozygous
|
||
sites indiscriminately alongside true duplications (see "Heterozygosity,
|
||
ploidy, and consensus-assembly inputs" below).</p>
|
||
<p><strong>Rejected: parsimony-based multiset pairing for multiplicity > 1.</strong> Rather
|
||
than dropping ambiguous loci, pair identical alleles between A and B first
|
||
(0-mutation explanation preferred), then take <code>min(unmatched_A, unmatched_B)</code>
|
||
as inferred SNP pairs. Rejected on two grounds: (1) circularity — selecting
|
||
pairs by minimal apparent divergence, then measuring divergence on those same
|
||
pairs, deflates the estimate by construction, not a neutral heuristic; (2)
|
||
the discriminating signal is a single base among 4 possible values, and the
|
||
flanks are <em>already</em> guaranteed identical for every candidate by
|
||
construction (that is how the locus was selected) — no information remains
|
||
in a k-mer window to tell which copy in B truly corresponds to which copy in
|
||
A once multiplicity > 1 on either side. Any pairing rule invents a
|
||
correspondence the data cannot support. Multiplicity > 1 is treated as
|
||
non-identifiable, not as a puzzle to solve with a heuristic.</p>
|
||
<h2 id="multi-genome-framing-family-as-pseudo-alignment-column">Multi-genome framing: family as pseudo-alignment column</h2>
|
||
<p><strong>Idea.</strong> Instead of resolving locus eligibility and correspondence one
|
||
genome pair at a time, treat a family as a column of a pseudo multiple
|
||
alignment across <em>all</em> genomes simultaneously: for each family, each genome
|
||
has either a net single-copy state (<code>A</code>/<code>C</code>/<code>G</code>/<code>T</code>, when the genome carries
|
||
exactly one of the 4 forms) or "missing" (<code>?</code>, multi-copy or absent). Flank
|
||
conservation (the <code>2m</code> bases fixed by construction) supplies positional
|
||
homology for free — the same role a real MSA would play, without alignment
|
||
software, gap penalties, or progressive-alignment approximations. Stacking
|
||
one such column per family, genomes as rows, produces a genuine SNP
|
||
pseudo-alignment matrix, not just a bag of pairwise distances.</p>
|
||
<p><strong>Precedent.</strong> This is the same principle behind reference-free
|
||
k-mer-based phylogenomics tools — SKA (Split K-mer Analysis, Harris 2018) and
|
||
kSNP: split the k-mer around a variable center, use flank identity to call
|
||
homologous columns across arbitrarily many genomes with no reference and no
|
||
MSA step, then feed the resulting pseudo-alignment to standard phylogenetic
|
||
tools. Landing on the same design independently is a good sign, not a
|
||
coincidence.</p>
|
||
<p><strong>Resolves the pairwise-correspondence problem, properly.</strong> The "Rejected:
|
||
parsimony-based multiset pairing" case above failed because, with only two
|
||
genomes' cardinalities to look at, there is no external constraint to justify
|
||
picking one correspondence between leftover alleles over another — <code>min(a,b)</code>
|
||
is a lower bound dressed up as a point estimate (see the follow-up discussion
|
||
on Felsenstein-style parsimony inconsistency: minimum-event explanations are
|
||
systematically biased low whenever homoplasy/multiplicity is real, not
|
||
noise-cancelling). With <code>N</code> genomes and many families jointly, the same
|
||
question can be answered the way real phylogenetics answers it: ancestral
|
||
state reconstruction / ML mapping over a tree estimated from the whole
|
||
column set. The tree supplies the missing constraint that two isolated
|
||
columns cannot — this is the principled way out, not a heuristic replacement
|
||
for one.</p>
|
||
<p><strong>Relation to what's already implemented.</strong> <code>KmerIndex::raw_snp_distance</code>
|
||
already computes, internally, per family, exactly this row — <code>single_form:
|
||
Vec<Option<u8>></code>, one entry per genome, <code>None</code> where ambiguous/absent —
|
||
before immediately collapsing it into pairwise <code>snp[i,j]</code>/<code>shared[i,j]</code>
|
||
tallies. The pivot this section proposes is small at the implementation
|
||
level: stop collapsing early, and surface the per-family row as a first-class
|
||
artifact (a <code>families x genomes</code> matrix). Pairwise raw p-distance becomes one
|
||
projection of that matrix (what's computed today), not the primary object;
|
||
downstream, the matrix itself could feed real phylogenetic tools (parsimony/
|
||
ML, e.g. RAxML/IQ-TREE-style) instead of only NJ/UPGMA on a homemade
|
||
pairwise-distance matrix.</p>
|
||
<p><strong>Caveat: column completeness shrinks with <code>N</code>.</strong> The probability that a
|
||
family's flanks stay intact simultaneously across all <code>N</code> genomes decays with
|
||
<code>N</code> (same ascertainment-bias mechanism as Bias 1 above, compounded over more
|
||
genomes) — fully-resolved columns (no <code>?</code> anywhere) become rare as more
|
||
genomes are added. Same missing-data situation any real multi-species
|
||
alignment faces, and phylogenetic tools already handle it well; the practical
|
||
implication is that columns should be allowed partial coverage (>=2 resolved
|
||
genomes, not unanimous) rather than requiring every genome to be net
|
||
single-copy at that locus.</p>
|
||
<h2 id="context-detectability-and-a-3-way-ordinal-distance-per-pair">Context, detectability, and a 3-way ordinal distance per pair</h2>
|
||
<p>Empirical follow-up to the pseudo-alignment idea above: <code>obikmer phylo
|
||
--snp</code> was run on a real 20-genome benchmark index and the resulting FASTA
|
||
fed to <code>raxml-ng</code>. Two problems surfaced, both traced back to conflating
|
||
distinct notions under one symbol.</p>
|
||
<p><strong>"Context", precisely.</strong> Sharing a central base between two genomes is not
|
||
just sharing a nucleotide — it is sharing a <strong>context</strong>: the <code>2m</code> flanking
|
||
bases, identical, which is a homology claim about that flanked window
|
||
(guaranteed non-coincidental by k-specificity, Bias 4 above), <em>not</em> a claim
|
||
about orthology or paralogy of the copy each genome carries. <code>A</code> opposite <code>C</code>
|
||
= same context, divergent centre. <code>A</code> opposite nothing = <strong>this context is
|
||
not observed in one of the two genomes</strong> — informative, not neutral.</p>
|
||
<p><strong>Why the IUPAC/DNA encoding used for the first <code>--snp</code> test was wrong.</strong>
|
||
Feeding IUPAC-coded ambiguity into a standard DNA model (<code>raxml-ng --model
|
||
GTR+G</code>) is a semantic mismatch: Felsenstein-pruning ML treats an ambiguous
|
||
tip as "exactly one true state, unknown which" (a uniform partial-likelihood
|
||
vector over compatible bases), not "these states are simultaneously
|
||
present". The two encodings look identical (same IUPAC letters) but the
|
||
software reads them backwards from what was intended — this invalidates the
|
||
literal branch lengths from that first experiment (topology-level groupings
|
||
by genus were still informative, see the worked example further down).</p>
|
||
<p><strong>Why <code>-</code> (absence) must not be scored as similarity, but also must not be
|
||
scored as a shared character between two absences.</strong> Two genomes both
|
||
lacking a context are not observed to resemble each other at that locus —
|
||
neither is observed to differ from the other either. It is a symmetric
|
||
non-observation, uninformative for that pair, and should contribute nothing
|
||
(not a small positive nor a small negative signal) to their distance. A
|
||
genome carrying a state (<code>A</code>) against one carrying none is a different case
|
||
entirely: informative, and should not be scored as neutral "missing data"
|
||
the way a generic DNA/ML pipeline would.</p>
|
||
<p><strong>Detectability vs existence — a deliberate simplification, accepted.</strong>
|
||
"Context not observed" conflates two different biological events: (1) true
|
||
loss of the locus, (2) the locus still exists but a mutation/indel <em>outside</em>
|
||
the centre, anywhere in the <code>2m</code> flanks, broke k-mer recognition. The design
|
||
adopts a rigorist stance on purpose: any flank-breaking mutation counts as
|
||
"this context no longer exists", full stop — because both causes (1) and (2)
|
||
independently require <em>at least</em> one more mutational event than a lone
|
||
central substitution would. This licenses treating "context absent in one of
|
||
the two genomes" as a <strong>lower bound</strong> on distance strictly greater than a
|
||
plain central SNP, without needing to know which of the two causes applies.
|
||
Coarser than a true event count, and accepted as such (fine substitution-type
|
||
modelling, e.g. transition/transversion weighting, is a secondary
|
||
refinement, not required for this to be useful).</p>
|
||
<p><strong>Resulting ordinal distance between two genomes at one family/context:</strong></p>
|
||
<table>
|
||
<thead>
|
||
<tr>
|
||
<th>Comparison</th>
|
||
<th>Distance</th>
|
||
<th>Meaning</th>
|
||
</tr>
|
||
</thead>
|
||
<tbody>
|
||
<tr>
|
||
<td>same centre (<code>A</code>/<code>A</code>)</td>
|
||
<td><code>0</code></td>
|
||
<td>identical</td>
|
||
</tr>
|
||
<tr>
|
||
<td>different centre, both single-copy (<code>A</code>/<code>C</code>)</td>
|
||
<td><code>1</code></td>
|
||
<td>plain central SNP</td>
|
||
</tr>
|
||
<tr>
|
||
<td>one genome has a state, the other has none</td>
|
||
<td><code>>1</code> (lower bound)</td>
|
||
<td>context undetectable in one genome — at least one extra mutational event, of unknown type</td>
|
||
</tr>
|
||
<tr>
|
||
<td>neither genome has any state (<code>∅</code>/<code>∅</code>)</td>
|
||
<td>excluded</td>
|
||
<td>symmetric non-observation, not comparable, contributes nothing</td>
|
||
</tr>
|
||
</tbody>
|
||
</table>
|
||
<p>This is a direct extension of <code>KmerIndex::raw_snp_distance</code> (<code>obikindex/src/siblings.rs</code>),
|
||
which today only implements the <code>0</code>/<code>1</code> rows and silently drops everything
|
||
else (including the informative <code>>1</code> row) rather than scoring it.</p>
|
||
<p><strong>Open, not yet resolved:</strong>
|
||
- Calibrating <code>>1</code> to a real number for tools expecting continuous distances
|
||
(NJ/UPGMA, ML branch lengths), rather than an arbitrary placeholder.
|
||
Natural route: estimate <code>p_hat</code> from the resolved (<code>0</code>/<code>1</code>) sites first,
|
||
then use the already-derived ascertainment formula (<code>P(usable window
|
||
showing a central SNP) = p * (1-p)^(2m)</code>, Bias 1 above) to derive a
|
||
model-consistent value for the <code>>1</code> bucket instead of guessing a constant.
|
||
- Where multi-copy/ambiguous states (the IUPAC case: a genome carrying more
|
||
than one form) fit into this ordinal scheme — plausibly also <code>>1</code> by the
|
||
same "at least one extra event" argument (a second form appearing is a
|
||
gain, itself an event), but not yet worked out.</p>
|
||
<p><strong>Practical alternative validated for the pseudo-alignment output itself</strong>
|
||
(orthogonal to the ordinal-distance question above, useful regardless of how
|
||
<code>>1</code> ends up calibrated): re-encode each family as 4 independent binary
|
||
presence/absence characters (<code>A</code>,<code>C</code>,<code>G</code>,<code>T</code> columns) instead of one IUPAC
|
||
column, feed to a <code>BIN</code>-type model instead of <code>DNA</code>. <code>∅</code> becomes an explicit
|
||
<code>0000</code> state (identity with another <code>0000</code>, not missing data) rather than a
|
||
gap — removes the semantic mismatch above by construction. Known cost,
|
||
accepted for now: a plain central substitution (<code>A</code> -> <code>C</code>) becomes 2 binary
|
||
flips (<code>1000</code> -> <code>0100</code>), overweighting substitutions relative to true
|
||
gain/loss events, and the 4 sub-characters of one family are not
|
||
statistically independent the way a generic <code>BIN</code> model assumes. A proper
|
||
fix (single 16-state alphabet, i.e. the powerset of <code>{A,C,G,T}</code>, with a
|
||
substitution-rate structure that respects the subset lattice rather than a
|
||
fully general 16x16 GTR-analogue) is very likely not expressible in
|
||
<code>raxml-ng</code>'s <code>MULTI</code> datatype as-is (Mk or fully-general rates only) and a
|
||
fully general 16-state rate matrix is almost certainly unidentifiable here
|
||
(states of cardinality >=3 are ~2% of sites in the benchmark run). Treated as
|
||
a longer-term research question, not a near-term implementation target.</p>
|
||
<h2 id="sankoff-parsimony-as-the-resolution-of-the-16-state-model-problem">Sankoff parsimony as the resolution of the 16-state model problem</h2>
|
||
<p>The "longer-term research question" just above (a 16-state alphabet — the
|
||
powerset of <code>{A,C,G,T}</code> — with a substitution structure that respects the
|
||
subset lattice) turns out to have a near-term answer, once the <em>unification</em>
|
||
question below is worked through.</p>
|
||
<p><strong>Distance methods (NJ/UPGMA/ME) vs. character methods (parsimony/ML): not a
|
||
deep philosophical divide, but a real practical distinction for this
|
||
project.</strong> Historically "phenetic" (characters -> distances -> tree) and
|
||
"cladistic" (characters -> tree directly) approaches were presented as
|
||
opposed schools; the modern view is a mathematical continuity, not a
|
||
dichotomy — Minimum Evolution (ME: find the tree minimizing total branch
|
||
length from a distance matrix) and Maximum Parsimony (MP: find the tree
|
||
minimizing total character-state changes) are both instances of "minimise a
|
||
global explanatory cost", and coincide under simple encodings (see Farris
|
||
1983, "The Logical Basis of Phylogenetic Analysis"; the MP/ME connection is
|
||
developed in the Minimum Evolution / Balanced Minimum Evolution literature,
|
||
e.g. Nei and colleagues — citations not independently re-verified here, flag
|
||
before quoting further). NJ's own agglomeration step already uses the whole
|
||
distance matrix jointly (the Q-matrix), not just the pair being merged — an
|
||
earlier claim in this discussion that distance methods are "blind" to
|
||
cross-taxon structure at every stage was too strong.</p>
|
||
<p>What <em>does</em> remain a real, structural distinction for this project: in a
|
||
character method, a given character's cost is <strong>re-evaluated per candidate
|
||
topology</strong> during tree search (the same family can cost 1 change under one
|
||
topology, 2 under another). In a pairwise-distance pipeline
|
||
(<code>raw_snp_distance</code> as it exists today), each family's contribution to
|
||
<code>d(i,j)</code> is computed <strong>once</strong>, independent of any candidate topology, before
|
||
NJ/UPGMA ever runs — so a question like "does this shared <code>∅</code> look like a
|
||
synapomorphy under topology T" can never be posed in that pipeline, for any
|
||
T. That question is only answerable by a method that tests candidate
|
||
topologies and re-scores characters under each — i.e. a character method.</p>
|
||
<p><strong>Sankoff parsimony directly resolves the <code>∅</code>/gain-loss/substitution
|
||
question, without the identifiability problem of a fitted 16-state model.</strong>
|
||
Sankoff's algorithm generalises Fitch parsimony to an arbitrary
|
||
user-supplied cost matrix between states (<code>obikseq</code>/<code>obikindex</code> would treat
|
||
each family as a <code>2^{4}</code>-state character, state = subset of <code>{A,C,G,T}</code>
|
||
observed in that genome, <code>∅</code> included as a real state, not a gap). The
|
||
previous 16-state idea failed specifically because <em>fitting</em> a full 16x16
|
||
rate matrix by ML is unidentifiable at this data volume; Sankoff sidesteps
|
||
that because the cost matrix is <strong>fixed a priori from domain knowledge</strong>, not
|
||
estimated — e.g. <code>c({A},{C}) = 1</code> (a substitution), <code>c({A},{A,C}) = 1</code> (a
|
||
gain), <code>c({A,C},{A}) = 1</code> (a loss), <code>c({A,C},{G,T}) = 2</code> (two changes) — no
|
||
estimation, no overparameterisation. This reframes gain/loss and central
|
||
substitution as two cost categories with independently chosen weights,
|
||
exactly the "two families of parameters" (<code>mu_substitution</code>, <code>mu_gain/loss</code>)
|
||
floated earlier in this discussion, now with an actual algorithmic home.</p>
|
||
<p><strong>Caveat, not blocking for this project's scope.</strong> Sankoff is still
|
||
parsimony: in principle exposed to Felsenstein's statistical-inconsistency
|
||
result under long-branch attraction (already invoked earlier against
|
||
<code>D_F = min(a,b)</code>) — parsimony and ML only provably coincide in the
|
||
short-branch regime. This is not a practical concern here because it is
|
||
exactly this estimator's declared target (closely related genomes, short
|
||
branches) — the regime where parsimony's known failure mode does not apply —
|
||
but worth stating explicitly as a scope guard rather than leaving it
|
||
implicit.</p>
|
||
<p><strong>Cheapest next experiment: don't write a Sankoff tree-search engine, use
|
||
one that exists.</strong> The hard part of a from-scratch implementation is not the
|
||
Sankoff DP itself (a straightforward dynamic program over a <em>fixed</em> tree) but
|
||
the topology search (SPR/NNI with incremental re-scoring) that comes for
|
||
free with <code>raxml-ng</code> on the ML side. <strong>TNT</strong> (Tree analysis using New
|
||
Technology, free, standard in morphological cladistics) already implements
|
||
Sankoff parsimony with a custom cost matrix plus topology search — the
|
||
family-state matrix (already close to what <code>--snp</code> produces, minus the
|
||
IUPAC/DNA-model mismatch) could be fed there directly, no new code required,
|
||
before considering a bespoke engine.</p>
|
||
<p><strong>The concrete comparison this unlocks:</strong> run both pipelines on the same
|
||
family data —
|
||
<code>k-mer families -> D_ij -> NJ/ME/UPGMA</code> (phenetic, what exists today) vs.
|
||
<code>k-mer families -> characters -> argmin_T Sankoff-cost(T)</code> (cladistic, via
|
||
TNT) — and compare the resulting topologies. Agreement would validate that
|
||
the pairwise-distance projection preserves the phylogenetic signal;
|
||
disagreement would pinpoint exactly what the projection to a single number
|
||
per pair loses. Not yet run.</p>
|
||
<h3 id="a-concrete-sankoff-cost-matrix-for-the-16-state-alphabet">A concrete Sankoff cost matrix for the 16-state alphabet</h3>
|
||
<p><strong>Why <code>c_gl</code> and <code>c_ctx</code> are the same constant, not two.</strong> An earlier version
|
||
of this design used two independent parameters: <code>c_gl</code> for an ordinary
|
||
gain/loss step between two nonempty states (e.g. <code>{A,C}->{A}</code>), and <code>c_ctx</code>
|
||
only for a state collapsing all the way to <code>∅</code>. That distinction doesn't
|
||
survive contact with what the index actually observes. Presence/absence
|
||
tracking never sees "the flanks" separately from "the centre" — it sees
|
||
whole, distinct, homologous 31-mers (two 31-mers sharing >=30 bases are
|
||
homologous by construction, though orthology vs. paralogy is undecidable
|
||
from that alone — settled earlier, see "Locus eligibility"). Losing the <code>A</code>
|
||
form of a family while <code>C</code> remains (<code>{A,C}->{C}</code>) and losing the last
|
||
remaining form (<code>{A}->∅</code>) are the <em>same kind of event</em>: a specific, complete
|
||
31-mer that used to be observed no longer is, because at least one of its
|
||
31 positions mutated. There is no separate "gain/loss of a still-recognised
|
||
allele" mechanism distinct from "loss of context" — both are exactly the
|
||
event <code>c_ctx(p)</code> (below) already computes the expected cost of. So: one
|
||
constant, <code>c_ctx</code>, used everywhere a member is gained or lost — <code>X -> ∅</code>
|
||
included, and (see next point) not even as a special case there.</p>
|
||
<p><strong>Set-edit-distance formula.</strong> For two states <code>X, Y ⊆ {A,C,G,T}</code>, split into
|
||
<code>seulement_X = X \ Y</code> (size <code>a</code>) and <code>seulement_Y = Y \ X</code> (size <code>b</code>).
|
||
Elements present in both cost nothing. Pair up to <code>min(a,b)</code> of the
|
||
remaining elements as <strong>substitutions</strong> (cheaper than treating them as an
|
||
unrelated loss + gain whenever <code>c_sub < 2*c_ctx</code>, which any sane parameter
|
||
choice satisfies); whatever is left over after pairing is a pure
|
||
<strong>gain/loss</strong>:</p>
|
||
<div class="highlight"><pre><span></span><code>cost(X,Y) = min(a,b)*c_sub + |a-b|*c_ctx
|
||
</code></pre></div>
|
||
<p>Worked examples (<code>c_sub = c_ctx = 1</code>): <code>{A}->{C}</code> = 1 (one substitution);
|
||
<code>{A}->{A,C}</code> = 1 (one gain, no substitution pair available since nothing is
|
||
only-in-Y that matches an only-in-X element after the shared <code>A</code> is
|
||
excluded); <code>{A,C}->{G,T}</code> = 2 (two substitution pairs, <code>A/C</code> vs <code>G/T</code>, both
|
||
same-size sets share nothing); <code>{A,C,G}->{A,C,T}</code> = 1 (<code>G</code>/<code>T</code> is the only
|
||
mismatched pair, <code>A,C</code> shared).</p>
|
||
<p><strong><code>∅</code> is <em>not</em> a flat-cost special case — a corrected position, reversing
|
||
an earlier draft of this design.</strong> An earlier version charged
|
||
<code>cost(X, ∅) = c_ctx</code> flat, independent of <code>|X|</code>, on the grounds that every
|
||
member of a family shares identical flanking sequence, so one mutation
|
||
breaking that context should unrecognise all of them at once rather than
|
||
<code>|X|</code> separate times. That argument doesn't survive comparison with how the
|
||
rest of the matrix already works: <code>{A,C,T} -> {A}</code> (losing two members,
|
||
one remaining) is charged <code>2*c_ctx</code> via the general formula above, with no
|
||
equivalent "maybe it was one shared event" discount — nothing distinguishes
|
||
that case from <code>{A,C} -> ∅</code> (losing two members, none remaining) other than
|
||
which state happens to be the target. Singling out <code>∅</code> for special
|
||
treatment was arbitrary, not principled: the uncertainty about "one event
|
||
or several" is identical in both cases, and the model already resolves it
|
||
uniformly elsewhere by simply counting elements. So <code>∅</code> is now an ordinary
|
||
node in the state graph below like any other, with no patch applied
|
||
afterwards:</p>
|
||
<div class="highlight"><pre><span></span><code>cost(X, ∅) = cost(∅, X) = |X| * c_ctx (via |X| single-element graph edges)
|
||
cost(∅, ∅) = 0
|
||
</code></pre></div>
|
||
<p><strong>Computing <code>c_ctx</code>.</strong> Not guessed — built from the value already derived
|
||
for the <code>>1</code> bucket in "Context, detectability, and a 3-way ordinal
|
||
distance per pair" above:</p>
|
||
<div class="highlight"><pre><span></span><code>c_ctx(p) = mean_sub_cost * [(2m*p) / (1 - (1-p)^(2m)) + p]
|
||
</code></pre></div>
|
||
<p><code>p</code> is <code>p_hat</code>, the calibrated per-site mutation probability (see
|
||
<code>calibrate_p_hat</code> / the "Experiment" section below); <code>m = (k-1)/2</code>, the
|
||
flank length on each side of the centre (<code>k=31</code> → <code>m=15</code>). Read the
|
||
bracketed term, well, term by term: <code>(1-p)^(2m)</code> is the probability that
|
||
<em>none</em> of the <code>2m</code> flanking positions mutated, so <code>1-(1-p)^(2m)</code> is the
|
||
probability that <em>at least one</em> did — i.e. the probability that this
|
||
specific 31-mer stops being observable at all, which is exactly the event
|
||
<code>c_ctx</code> prices. <code>2m*p</code> is the unconditional expected mutation count over
|
||
those <code>2m</code> positions. Their ratio is the conditional expectation
|
||
<code>E[mutations | at least one occurred]</code> — provable directly: for any
|
||
nonnegative integer random variable <code>X</code>, <code>X · 1{X>=1} = X</code> always (both
|
||
sides are <code>0</code> when <code>X=0</code>, both are <code>X</code> otherwise), so
|
||
<code>E[X | X>=1] = E[X·1{X>=1}]/P(X>=1) = E[X]/P(X>=1)</code>, which is
|
||
<code>(2m*p)/(1-(1-p)^(2m))</code> here. The trailing <code>+p</code> adds the (much smaller,
|
||
first-order, not itself a conditional expectation) marginal contribution
|
||
of the centre position's own mutation probability.</p>
|
||
<p><strong>The <code>mean_sub_cost</code> multiplier — a correction, not part of the original
|
||
derivation.</strong> The bracketed term is a <em>count</em> of expected mutations, not a
|
||
cost — an earlier version used it directly as <code>c_ctx</code>, implicitly pricing
|
||
every one of those mutations at a flat <code>1</code> regardless of type. That stopped
|
||
being defensible once substitution costs were calibrated per base-pair
|
||
category (<code>substitution_costs_from_tally</code>): transitions are markedly more
|
||
frequent than transversions in this project's own data (roughly 5-6x, e.g.
|
||
130k vs. 20-28k observed instances per category in the eubacteria run
|
||
below) and correspondingly cheaper, so "one mutation" isn't worth a flat
|
||
unit. <code>mean_sub_cost</code> (<code>mean_substitution_cost</code>) is the empirical average
|
||
substitution cost, weighted by each category's observed frequency — it
|
||
converts the bracketed term's expected mutation <em>count</em> into an actual
|
||
expected <em>cost</em>. Concretely, in the eubacteria run: the bracketed term
|
||
alone is <code>~1.17</code> (close to <code>1</code>, i.e. context loss is <em>usually</em> attributable
|
||
to a single mutation, since <code>p_hat</code> is small); <code>mean_sub_cost ~= 1.48</code>
|
||
(pulled up from the cheapest transition cost of <code>1.0</code> by the substantial
|
||
minority of transversions in the mix); final <code>c_ctx ~= 1.73</code> — sensibly
|
||
between the cheapest transition (<code>~1.0</code>) and the transversion costs
|
||
(<code>~2.5</code>-<code>2.9</code>), rather than coinciding almost exactly with the cheapest
|
||
transition purely by construction accident, which is what the unmultiplied
|
||
version did and which is what exposed this gap in the first place.</p>
|
||
<p><strong>Better construction method: shortest path in a small state graph, not the
|
||
closed-form formula directly.</strong> Build a graph on all 16 states, <code>∅</code>
|
||
included, with two edge types — substitution edges between same-cardinality
|
||
sets differing by one element (weight <code>c_sub</code>, or <code>c_ts</code>/<code>c_tv</code> if split
|
||
further below), and gain/loss edges between sets whose cardinality differs
|
||
by one (weight <code>c_ctx</code>, per the unification above; <code>∅</code> connects to each
|
||
singleton state this way, being a strict subset of it) — then define
|
||
<code>cost(X,Y)</code> as shortest-path distance in that graph, uniformly for every
|
||
<code>X,Y</code> including <code>∅</code>, precomputed once (16 nodes, trivial) into a dense
|
||
16x16 matrix before feeding it to Sankoff/TNT. Verified equivalent to the
|
||
closed-form formula above in the
|
||
uniform-cost case (checked by hand on <code>{A}->{C,G,T}</code>: both give
|
||
<code>c_sub + 2*c_ctx</code>). The graph construction is not just a reformulation for
|
||
its own sake: it is the version that generalises correctly once substitution
|
||
costs stop being uniform (next point) — the closed-form's <code>min(a,b)</code>
|
||
counting silently assumes <em>any</em> pairing costs the same, which breaks the
|
||
moment <code>c_sub</code> depends on which two bases are involved.</p>
|
||
<p><strong>Substitution refinement — implemented as the full 6-category symmetric
|
||
matrix, not just Ts/Tv.</strong> <code>c_sub</code> was originally going to split into just
|
||
<code>c_ts</code> (A<->G or C<->T) and <code>c_tv</code> (the other four pairs) — already
|
||
well-defined in canonical space (see "Canonical invariance" above: a
|
||
transition maps to a transition, a transversion to a transversion,
|
||
regardless of orientation). A symmetric cost matrix allows finer
|
||
resolution than that 2-category split, though: each of the 6 distinct
|
||
undirected base pairs (AC, AG, AT, CG, CT, GT) can be costed
|
||
independently, with no grouping at all — <code>SankoffWeights::sub_cost</code> is a
|
||
full <code>[[f64;4];4]</code> table (<code>SankoffWeights::ts_tv</code> remains as a convenience
|
||
constructor for the coarser 2-category case, if ever wanted). This turns
|
||
"pick <code>min(a,b)</code> substitution pairs" into a genuine (tiny, <=4 elements per
|
||
side, trivially enumerable) minimum-cost bipartite matching problem instead
|
||
of a plain count — the state-graph shortest-path construction handles this
|
||
automatically, no separate logic needed. Calibration is not a guess
|
||
either: <code>KmerIndex::base_pair_tally</code> collects the pooled <code>(centre_i,
|
||
centre_j)</code> distribution over resolved (SNP) sites — restricted to the same
|
||
ratio-ceiling-included genome pairs as <code>p_hat</code>'s own calibration, for the
|
||
same saturation-exclusion reason — and <code>substitution_costs_from_tally</code>
|
||
turns that into <code>cost(a,b) = -ln(rate(a,b))</code>, normalised so the most
|
||
frequent category costs <code>1.0</code> (the standard generalised-parsimony
|
||
step-weighting heuristic, generalised from 2 categories to 6). A second
|
||
full pass over the annex is required for this (see <code>base_pair_tally</code>'s own
|
||
docs for why it can't share <code>raw_snp_distance</code>'s single pass) — both now
|
||
share one traversal helper, <code>scan_family_pairs</code>, rather than duplicating
|
||
the per-family reconciliation logic.</p>
|
||
<p><strong>Caveat carried over from the <code>D_F = min(a,b)</code> rejection earlier:</strong> this
|
||
cost matrix is only valid as <strong>Sankoff step-cost input</strong>, re-evaluated for
|
||
every branch of every candidate topology during tree search. Reusing
|
||
<code>cost(leaf_A, leaf_B)</code> directly as a standalone pairwise distance (bypassing
|
||
the tree) would reintroduce the exact circularity already rejected — the
|
||
<code>min(a,b)</code> pairing here is a locally-defined edit distance between two
|
||
states, not a claim about the true evolutionary history between two
|
||
specific genomes.</p>
|
||
<p><strong>Feasibility confirmed:</strong> TNT's <code>costs</code> command accepts custom step
|
||
matrices for multistate characters, so this whole construction (16x16
|
||
matrix derived from the state graph, <code>c_ts</code>/<code>c_tv</code>/<code>c_ctx</code> as the three
|
||
tunable parameters) is directly usable there — no new tooling required
|
||
before testing it.</p>
|
||
<h3 id="experiment-tnt-run-on-real-data-2026-08-11">Experiment: TNT run on real data (2026-08-11)</h3>
|
||
<p><strong>Status:</strong> validated at genus/family/order scale within Bacteria; not
|
||
informative across domains with this character type. Exploratory only —
|
||
run entirely outside the repo (<code>/tmp/tnt_run</code>, TNT installed locally under
|
||
<code>TNT/</code>), no new Rust code. Kept here as the record of what was learned.</p>
|
||
<p><strong>Pipeline.</strong> <code>obikmer phylo --snp</code> emits one IUPAC-coded pseudo-alignment
|
||
row per genome (<code>snp_pseudo_alignment</code>, one column per family with
|
||
<code>family_size() >= 2</code>). A small external Python script decodes IUPAC back to
|
||
the 16-state bitmask, applies the set-edit-distance formula above, and
|
||
emits a complete TNT script (<code>xread</code> matrix + <code>smatrix</code> step-matrix + <code>hold</code>/
|
||
<code>mult</code> search). No new Rust code was needed for this pass.</p>
|
||
<p><strong>Calibration ("quick option").</strong> Rather than modifying
|
||
<code>write_raw_snp_distance_csv</code> to emit raw counts, <code>c_ctx</code> was approximated as
|
||
the unweighted mean of the pairwise <code>snp/(snp+shared)</code> ratios already
|
||
present in <code>rawsnp.csv</code>, restricted to the relevant taxon subset. Used
|
||
values: <code>mu = gamma = 10</code>, <code>c_ctx = 18</code> (this run predates the <code>c_gl</code>/<code>c_ctx</code>
|
||
unification above — <code>gamma</code> here is what <code>c_gl</code> was called before it turned
|
||
out to just be <code>c_ctx</code>; ratio <code>c_ctx/mu = 1.8</code>, consistent
|
||
with the observed pairwise ratios for genomes within Enterobacteriaceae/
|
||
Eubacteria, which cluster around 1.5-2 and barely move as the taxon set
|
||
widens — the context signal is stable, not sensitive to which subset is
|
||
chosen). The more principled route (raw counts, weighted <code>p_hat</code>, Ts/Tv
|
||
split) remains a follow-up, not yet done.</p>
|
||
<p><strong>Run 1 — Enterobacteriaceae (11 taxa: 4 <em>E. coli</em>, 4 <em>Salmonella enterica</em>,
|
||
3 <em>Klebsiella pneumoniae</em>).</strong> All three genera recovered as monophyletic.
|
||
Initial read of the exported (unrooted, TNT/Nexus <code>[&U]</code>) tree as showing a
|
||
genus-arrangement disagreement with known systematics (<em>Escherichieae</em>:
|
||
<em>Escherichia</em>+<em>Salmonella</em> sister vs. more distant <em>Klebsielleae</em>) was
|
||
<strong>wrong</strong> — diagnosed via <code>force = (taxa);</code> monophyly-constraint test
|
||
(identical score constrained vs. free ⟹ no real disagreement). Root cause of
|
||
the misreading: with exactly 3 clades and no outgroup, an unrooted tree has
|
||
only <strong>one possible topology</strong> (a single trifurcation) — there is no
|
||
internal arrangement to get right or wrong. This run cannot test
|
||
inter-genus relationships at all; it can only test intra-genus monophyly
|
||
(which held).</p>
|
||
<p><strong>Run 2 — Eubacteria (18 taxa: run 1 + <em>Acidobacterium capsulatum</em>,
|
||
<em>Opitutus terrae</em>, <em>Bacillus subtilis</em>, <em>Shouchella clausii</em>, <em>Wolbachia</em>
|
||
endosymbiont, <em>Proteus mirabilis</em>, <em>Yersinia ruckeri</em>).</strong> The
|
||
Enterobacteriaceae substructure from run 1 is reproduced identically, now
|
||
correctly rooted by real outgroups, resolving the tribal arrangement left
|
||
undetermined in run 1: <em>Escherichia</em>+<em>Salmonella</em> sister, <em>Klebsiella</em> more
|
||
distant — matching known systematics. Outgroup placement: <code>(Bacillus,
|
||
Shouchella)</code> sister pair (Firmicutes/<em>Bacillales</em>) splits from all
|
||
Proteobacteria — a correct phylum-level split; <code>Wolbachia</code>
|
||
(Alphaproteobacteria) splits from the Gammaproteobacteria block
|
||
(<code>Proteus</code>, <code>Yersinia</code>, Enterobacteriaceae) — a correct class-level split.
|
||
Lower confidence: the fine nested order <code>(Proteus, (Yersinia,
|
||
Enterobacteriaceae))</code> and the relative position of <code>Acidobacterium</code> vs.
|
||
<code>Opitutus</code> — plausible, not independently verified against current
|
||
Enterobacterales family-level literature.</p>
|
||
<p><strong>Run 3 — full domain set (20 taxa: run 2 + <em>Candidozyma auris</em> [yeast] and
|
||
<em>Saccharolobus islandicus</em> [archaeon]).</strong> The bacterial clade from run 2 is
|
||
reproduced <strong>unchanged and intact</strong> — a real robustness signal, the method
|
||
does not fragment the ingroup when unrelated deep taxa are added. But the
|
||
result carries <strong>no information on Bacteria/Archaea/Eukarya relationships</strong>:
|
||
with exactly one eukaryote, one archaeon and one bacterial clade, the
|
||
unrooted tree is again forced into the single 3-clade trifurcation from run
|
||
1's caveat — there is no second representative of either outgroup domain to
|
||
resolve internal arrangement, so nothing about their relative position can
|
||
be read from the topology (a "ladder" ordering in the exported tree is
|
||
serialization, not signal). Independently, <code>rawsnp.csv</code> shows <em>why</em> this
|
||
character type cannot reach further: pairwise ratios involving
|
||
<em>Candidozyma</em>/<em>Saccharolobus</em> are almost all <code>NA</code> (no central-position
|
||
family shared at all) or saturated at <code>1.0</code> (every shared family differs) —
|
||
central-position families require literal 31 bp context conservation, which
|
||
simply does not survive domain-level divergence. Cross-domain placement
|
||
would need conserved-marker characters (rRNA, ribosomal proteins), not this
|
||
estimator.</p>
|
||
<h3 id="native-sankoff-tnt-phyg-export-2026-08-12-superseding-the-external-scripts-above">Native <code>--sankoff --tnt</code>/<code>--phyg</code> export (2026-08-12), superseding the external scripts above</h3>
|
||
<p>The ad hoc Python glue from the previous section is superseded: <code>obikmer
|
||
distance --sankoff</code> now calibrates the matrix natively (<code>p_hat</code>, 6-category
|
||
substitution costs, <code>c_ctx</code> weighted by <code>mean_sub_cost</code> — see the worked
|
||
example above) and <code>--tnt</code>/<code>--phyg</code> each write a ready-to-run script from
|
||
it, no external script needed. <code>∅</code> is an ordinary 16th state throughout
|
||
(never <code>-</code>), specifically to avoid gap-semantics confusion in downstream
|
||
tools — see "A concrete Sankoff cost matrix" above for why.</p>
|
||
<p><strong>TNT (<code>--tnt</code>).</strong> <code>write_sankoff_tnt</code> (<code>obikmer/src/cmd/phylo/mod.rs</code>)
|
||
recodes to TNT's own <code>0-9A-F</code> xread alphabet (its default reader rejects
|
||
the wider IUPAC set otherwise), scales and rounds costs to integers
|
||
(<code>smatrix</code>/<code>cost</code> reject decimals), then re-runs integer Floyd-Warshall on
|
||
the rounded matrix (<code>scaled_metric_matrix</code>) — independently rounding each
|
||
cell of an already-metric real-valued matrix can break the triangle
|
||
inequality (e.g. two real costs of <code>1.734</code> round to <code>173</code> each, summing to
|
||
<code>346</code>, while their real sum <code>3.468</code> rounds to <code>347</code>), which TNT otherwise
|
||
silently "fixes" itself with an unreproducible correction. Verified against
|
||
the real 20-genome benchmark index: zero triangle-inequality violations
|
||
after the fix, TNT loads the file without its "triangle inequality
|
||
violated... Fixed" warning.</p>
|
||
<p>Two syntax facts worth recording because they're wrong in intuitive
|
||
guesses and contradicted actual TNT behavior when tested: TNT's plain
|
||
command stream has <strong>no comment syntax</strong> of its own — <code>/* */</code> and <code>[ ]</code>
|
||
only work inside the (separately-enabled, off by default) macro scripting
|
||
language, and error with "No command!" otherwise. The working substitute
|
||
is <code>quote TEXT ;</code> (prints the text, doesn't affect parsing) — but the text
|
||
itself can't contain a literal <code>;</code> (TNT's universal terminator); the
|
||
manual's own escape (<code>.,</code>) exists but the script here just avoids
|
||
semicolons in the text instead.</p>
|
||
<p>The default search command embedded in the script is <code>mult</code> (traditional:
|
||
random addition sequences + TBR), not <code>xmult</code> (New Technology search:
|
||
ratchet/drift/tree-fusion). <code>xmult</code> with TNT's default <code>mxram</code> (16 MB,
|
||
must be set <em>before</em> <code>xread</code> if changed) ran out of RAM on the real
|
||
908k-character dataset ("<code>xmult - out of ram</code>"); <code>mult</code> does not, matching
|
||
what had already been validated by hand outside this session.</p>
|
||
<p><strong>PhyG (<code>--phyg</code>).</strong> <code>write_sankoff_phyg</code> writes a <code>tcm:</code> custom-alphabet
|
||
matrix (same scale+round+metric-closure treatment as TNT) and reuses
|
||
<code>--sankoff</code>'s own <code>_sankoff.fasta</code> as-is via <code>prefasta:</code> — PhyG's <code>tcm:</code>
|
||
alphabet is read from the matrix file's own first line, so (unlike TNT) no
|
||
recoding is needed. PhyG auto-adds its own indel/gap state as an
|
||
<code>(n+1)</code>-th row/column of the tcm; inert here since the alignment encodes
|
||
absence as <code>0</code>, never <code>-</code>.</p>
|
||
<p><code>report("file", newick, overwrite)</code> — exactly as shown in PhyG's own
|
||
manual — triggers <code>Unrecognized/missing report option ... defaulting to
|
||
'graphs'</code> on the locally installed binary (1.3, commit <code>3c1a1fa</code>); the
|
||
working form adds <code>graphs</code> explicitly: <code>report("file", graphs, newick,
|
||
overwrite)</code>. Manual/binary mismatches like this (also true of <code>criterion:</code>
|
||
— the binary accepts <code>parsimony</code>/<code>ml</code>/<code>pmdl</code>, the manual instead documents
|
||
<code>mapa</code>/<code>ncm</code>/<code>parsimony</code>/<code>pmdl</code>/<code>si</code>) mean command syntax against this PhyG
|
||
build should be verified empirically, not trusted from the PDF alone.
|
||
<code>instances:N</code> (not a separate CPU flag) is what actually parallelises the
|
||
search across cores — PhyG uses all physical cores by default but only
|
||
across as many instances as are running, so raise it to the physical core
|
||
count to use them all (the CLI-level <code>+RTS -NX -RTS</code> flag also exists but
|
||
controls something else: capping/limiting cores, not requesting more).</p>
|
||
<p>Both scripts share one <code>--sankoff-cost-scale</code> (default <code>100</code>), not two
|
||
separate flags — they scale the same calibrated matrix for the same
|
||
reason (integer-only cost commands) and no PhyG-specific
|
||
accumulator-width constraint was ever found to justify a different
|
||
default from TNT's (TNT: hinted 32-bit accumulators in its own manual;
|
||
PhyG: no such hint found — Haskell's native <code>Int</code> is typically 64-bit).</p>
|
||
<p><strong>Open problem: PhyG reports all branch lengths as <code>0.0</code>.</strong> The graph-level
|
||
parsimony cost is correct (<code>3.3286×10⁸</code> on the real dataset, consistent
|
||
with TNT's <code>328574911</code> on the same calibrated matrix), but every individual
|
||
edge in the exported Newick shows <code>:0.0</code>, with the total cost only ever
|
||
shown as a whole-tree annotation (<code>[3.32860377e8]</code>). Not fixed, not fully
|
||
diagnosed — PhyG's manual describes per-edge branch length as computed by
|
||
ancestral-state (HTU) backtracking, well documented for sequence/standard
|
||
character types, but nothing found (the term "Sankoff" doesn't even appear
|
||
in the manual) confirming this backtracking is wired up for a custom
|
||
<code>tcm:</code> matrix character. Switching <code>criterion:</code> to a likelihood-family
|
||
option (<code>ml</code>, or the manual's <code>mapa</code>/<code>ncm</code>/<code>si</code>) was considered as a
|
||
possible fix but is very unlikely to be one: those criteria are
|
||
information-theoretic reparametrisations of the <em>same</em> step-counting
|
||
machinery as parsimony (<code>ncm</code> in particular is known in the literature to
|
||
be numerically equivalent to weighted parsimony), not classical
|
||
continuous-time-Markov ML with a real rate matrix — so they wouldn't
|
||
change how branch length is attributed per edge either.</p>
|
||
<p><strong>Export format note (not a bug in the generator).</strong> Neither script's
|
||
<code>.tre</code> output opens in PearTree (FigTree's successor) via File > Open —
|
||
association/Launch-Services quirks were ruled out (the file was opened
|
||
directly through the app, not by double-click). Likely cause, not yet
|
||
confirmed: TNT's export is a minimal NEXUS <code>begin trees;</code> block with no
|
||
preceding <code>Taxa</code> block and bare numeric (untranslated) leaf labels; PhyG's
|
||
is multiple raw Newick trees concatenated with no NEXUS wrapper at all
|
||
plus a trailing <code>[cost]</code> bracket tag after the root label. Both differ
|
||
from a "normal" single-tree, fully-declared NEXUS file; this is PhyG/TNT's
|
||
own export format, not something <code>write_sankoff_tnt</code>/<code>write_sankoff_phyg</code>
|
||
could fix without post-processing the <em>other</em> program's output after the
|
||
fact.</p>
|
||
<h3 id="next-direction-genuine-ml-branch-lengths-not-parsimony-open-2026-08-12">Next direction: genuine ML branch lengths, not parsimony (open, 2026-08-12)</h3>
|
||
<p>Decided: parsimony (the whole <code>--sankoff</code>/<code>--tnt</code>/<code>--phyg</code> pipeline above)
|
||
is a stopgap, not the destination. The goal is maximum likelihood with
|
||
real, calibrated branch lengths (expected substitutions/site), which
|
||
parsimony step-counts were never going to give directly (see the open
|
||
"branch lengths are <code>0.0</code>" problem above — even if fixed, TNT/PhyG-style
|
||
parsimony branch length is a step count, not a continuous ML estimate).</p>
|
||
<p><strong>Model choices, settled:</strong>
|
||
- <strong>The exchangeability <code>R</code> is symmetric; the rate matrix <code>Q</code> is not.</strong>
|
||
(Superseded an earlier, wrong framing here that treated "symmetric
|
||
model" as one thing — see the resolution below on <code>R</code> vs <code>Q</code> vs <code>π</code> for
|
||
the full reasoning.) <code>R(a,b) = R(b,a)</code> because <code>BasePairTally</code> never
|
||
captured direction — a fact about the data, not a modelling choice.
|
||
<code>Q(a,b) = R(a,b)·π_b</code> is asymmetric whenever the real state frequencies
|
||
<code>π</code> are (which they are, empirically) — biology drives this via <code>π</code>,
|
||
not via <code>R</code>.
|
||
- <strong><code>∅</code> stays an ordinary 16th state</strong>, as already established for
|
||
TNT/PhyG — same reasoning applies to any ML tool: encode as a real
|
||
alphabet symbol, never as <code>-</code>/gap, or the RAxML-era failure (empty set
|
||
silently treated as missing data) repeats.</p>
|
||
<p><strong>Stationary frequencies for the 16 states — resolved (2026-08-12).</strong>
|
||
A CTMC needs a rate matrix <code>Q</code>, generally asymmetric. <code>Q(i,j) = R(i,j) ·
|
||
π_j</code>, where <code>R</code> (exchangeability) is symmetric and <code>π</code> (stationary
|
||
frequencies) need not be — this factoring is what makes <code>Q</code> reversible
|
||
(satisfies detailed balance, <code>π_i·Q(i,j) = π_j·Q(j,i)</code>) for <em>any</em> <code>π</code>, not
|
||
just uniform, as long as <code>R</code> is symmetric. Two separate, both-easy
|
||
quantities, not one hard inverse problem:
|
||
- <strong><code>R</code> is already calibrated</strong>: <code>sub_cost</code> (<code>-ln(observed rate)</code> per pair,
|
||
from <code>BasePairTally</code>) <em>is</em> <code>R</code> up to a log transform — recover it as
|
||
<code>R(a,b) = exp(-sub_cost(a,b))</code>. Symmetric by construction, because the
|
||
tally itself never distinguished direction (unordered-pair counts only)
|
||
— not a modelling choice, a fact about what the data can say.
|
||
- <strong><code>π</code> is a direct count</strong>: empirical marginal frequency of each of the
|
||
16 states across the whole alignment (same kind of scan already used to
|
||
confirm <code>N</code> occurs 1383 times in the real 20-genome benchmark). With
|
||
~908k sites × 20 genomes, the counts are large enough that this is
|
||
precise on its own — no need to spend ML degrees of freedom
|
||
re-estimating it via IQ-TREE's <code>+FO</code>.</p>
|
||
<p>Checked and ruled out along the way: IQ-TREE's <code>+F</code> (empirical, "compute
|
||
from the alignment") does <strong>not</strong> work as a shortcut for this — for a
|
||
custom-file morphology model, <code>readParameters</code> always requires the file's
|
||
own frequency line unconditionally; omitting it and passing <code>+F</code> instead
|
||
just fails (<code>ERROR: State frequencies could not be read</code>). <code>π</code> has to be
|
||
computed by <code>obikmer</code> and written into the file, not left to IQ-TREE.</p>
|
||
<p>Net effect: no free-rate ML estimation needed for this piece at all (the
|
||
mistaken assumption that motivated most of this discussion — that
|
||
building an asymmetric <code>Q</code> a priori would require solving a linear system
|
||
from <code>Q</code> itself — doesn't apply, because <code>R</code>, the only piece that's
|
||
genuinely hard to get directionally, is symmetric and already in hand).</p>
|
||
<p><strong>Candidate tool: IQ-TREE</strong>, because it supports user-defined multistate
|
||
models (unlike RAxML's <code>MULTI</code> data type, which is limited to the
|
||
equal-rate Mk model and can't take a custom rate matrix at all — a genuine
|
||
tool limitation, not a gap-symbol encoding problem this time). IQ-TREE 3
|
||
(3.0.1) is now installed locally (Homebrew, <code>iqtree3</code>).</p>
|
||
<p><strong>IQ-TREE custom-model format — verified empirically against the local
|
||
binary (2026-08-12).</strong> The web docs' <code>-mdef</code> NEXUS <code>begin models;
|
||
frequency NAME = ...; model NAME = ...; end;</code> mechanism (initially assumed
|
||
to apply directly, see history below) turned out to be for <strong>named
|
||
components used inside <code>MIX{...}</code>/<code>FMIX{...}</code> mixture models only</strong> — it
|
||
does <strong>not</strong> apply to a single, non-mixture custom morphology matrix, and
|
||
using it that way fails (<code>ERROR: File not found <name></code> — traced in
|
||
IQ-TREE 3's own source, <code>model/modelmorphology.cpp</code>: any <code>-m</code> string that
|
||
isn't <code>MK</code>/<code>ORDERED</code>/<code>GTR</code>/<code>GTRX</code> is passed straight to
|
||
<code>ModelMarkov::readParameters()</code>, which opens it <strong>as a literal file path</strong>,
|
||
never consulting the <code>-mdef</code> models block at all for this data type).</p>
|
||
<p><strong>The confirmed working recipe</strong> (built a tiny 5-taxon/3-state toy dataset
|
||
and rate file, ran it end to end with <code>iqtree3</code>, got a real ML tree with
|
||
non-zero branch lengths and an optimized log-likelihood — ground truth,
|
||
not documentation):
|
||
- No <code>-mdef</code> needed. Write one plain file (any name) containing, as
|
||
whitespace/newline-separated numbers, in order: the <strong>lower-triangular
|
||
rate matrix</strong> (<code>N(N-1)/2</code> values, PAML row-major order — for 16 states,
|
||
120 values, the same count and layout already produced for TNT's
|
||
<code>smatrix</code>), immediately followed by the <strong>N state frequencies</strong> on the
|
||
same stream (no header, no blank line required — confirmed by reading
|
||
<code>ModelMorphology::readRates</code>/<code>ModelMarkov::readStateFreq</code> directly, which
|
||
just pull tokens off the stream in sequence).
|
||
- Invoke with <code>-m <path-to-that-file>+ASC</code> (<code>+ASC</code> for the no-constant-site
|
||
correction, as before). An explicit <code>+F{f1,...,fN}</code> on the command line
|
||
overrides the file's own frequency line if given (confirmed in
|
||
<code>ModelMorphology::init</code>) — useful once real calibrated stationary
|
||
frequencies exist, a placeholder equal-frequency line works meanwhile
|
||
(the still-open gap noted above).
|
||
- <code>--seqtype MORPH</code> (alphabet <code>0</code>-<code>9</code>,<code>A</code>-<code>V</code>, ≤31 states) — reuse the same
|
||
<code>0-9A-F</code> recoding already built for TNT (<code>TNT_STATE_SYMBOL</code>).</p>
|
||
<p><strong>Risk, confirmed, precisely characterised, and resolved by design
|
||
(2026-08-12).</strong>
|
||
<code>--seqtype MORPH{16}</code> <strong>does not force the state count</strong> for real ML
|
||
analysis — tested directly (<code>--seqtype MORPH{4}</code> on the 3-symbol toy
|
||
alignment gave the byte-for-byte identical 3-state result as no <code>{4}</code> at
|
||
all) and confirmed in source: the value it sets
|
||
(<code>params.alisim_num_states_morph</code>, <code>utils/tools.cpp</code>) is consumed only by
|
||
the <code>--alisim</code> simulator; the main analysis path always calls
|
||
<code>getDataBlockMorphStates</code>/an equivalent scan (<code>alignment/alignment.cpp</code>),
|
||
for both FASTA/PHYLIP and NEXUS input (a NEXUS <code>symbols=</code> declaration
|
||
doesn't change this either — checked, same code path). No CLI flag or
|
||
NEXUS declaration overrides it.</p>
|
||
<p>The precise rule (from <code>getDataBlockMorphStates</code>, <code>alignment.cpp:1058</code>):
|
||
<code>N</code> = <strong>one plus the highest state ordinal actually observed anywhere in
|
||
the alignment</strong>, ordinal being the symbol's position in IQ-TREE's own
|
||
fixed table <code>"0123456789ABCDEFGHIJKLMNOPQRSTUV"</code> — not a count of distinct
|
||
symbols seen. So the risk is narrower than "any missing symbol breaks it":
|
||
concretely, it's whether the symbol mapped to state index 15 (<code>F</code> in the
|
||
<code>0-9A-F</code> recoding already used for TNT, i.e. our <code>N</code> = "all four bases
|
||
ambiguous") occurs <strong>at least once anywhere</strong> in the real alignment — if
|
||
it does, <code>N</code> correctly comes out to 16 regardless of which lower-index
|
||
symbols (including <code>0</code>/<code>∅</code>) are rare or absent; if it doesn't, <code>N</code> silently
|
||
undercounts and misaligns every value in a 16-entry rate/frequency file,
|
||
with no error to catch it. <strong>Checked against the real biological alignment</strong>
|
||
(<code>/tmp/msg_test/eub_sankoff.fasta</code>, the 20-genome benchmark index): the
|
||
symbol <code>N</code> (IUPAC "all four bases ambiguous," state index 15) occurs 1383
|
||
times across 13 of the 20 sequences — present, so this specific real
|
||
dataset is not at risk. Still worth a real, general presence check inside
|
||
<code>obikmer</code> before this is wired in, rather than assuming every future
|
||
dataset will have <code>N</code> too (nothing in IQ-TREE would catch it if not).</p>
|
||
<p><strong>Resolution: subset + compact-renumber, not rely on all 16 appearing.</strong>
|
||
Since IQ-TREE always infers <code>N</code> from the alignment's own content and
|
||
nothing overrides that, the fix is to make the file <code>obikmer</code> writes match
|
||
that inference <em>by construction</em>, for every run, rather than hope the
|
||
16th (or any particular) state happens to occur:
|
||
1. Scan the real alignment for which of the 16 canonical states actually
|
||
occur anywhere (not per-column — anywhere in the whole alignment).
|
||
2. Renumber the occurring states to a <strong>compact, consecutive</strong> <code>0..k-1</code>
|
||
range, preserving their relative order (the original bitmask/<code>STATE_
|
||
SYMBOL</code> ordering) — not just filtering, since a <em>gap</em> in the ordinal
|
||
sequence (e.g. keeping states <code>{0,1,2,4}</code> numbered as-is instead of
|
||
<code>{0,1,2,3}</code>) reproduces the exact same "highest observed ordinal"
|
||
miscount this was meant to fix.
|
||
3. Recode the alignment itself with this new compact <code>k</code>-symbol alphabet
|
||
(same recoding mechanism already used for TNT's <code>0-9A-F</code>, just over a
|
||
possibly-smaller symbol set).
|
||
4. Extract the matching <code>k×k</code> submatrix (rows/columns for the kept
|
||
states only) from the full calibrated 16×16 cost matrix, in the same
|
||
lower-triangular order the rate-matrix file needs — and, later, the
|
||
matching <code>k</code>-length subset of stationary frequencies once those are
|
||
calibrated (still the open gap noted earlier in this section).</p>
|
||
<p>Consequence, and why nothing is lost: a state that never occurs in a given
|
||
alignment can, by definition, never contribute a transition to score in
|
||
that same alignment — dropping it from that run's matrix costs nothing.
|
||
The subset (and therefore <code>k</code>) can differ from one dataset/run to the
|
||
next; this has to be done freshly per alignment, not computed once and
|
||
reused.</p>
|
||
<p><em>(Superseded reasoning, kept for the record: the <code>-mdef</code> NEXUS route
|
||
below was the original plan, based on IQ-TREE's own web documentation for
|
||
protein mixture models, before the local install allowed testing it —
|
||
<code>GTRX</code> combined with a <code>-mdef</code>-referenced custom model, <code>+Fname</code> frequency
|
||
reference. Both pieces exist and parse without error individually, but
|
||
<code>GTRX</code>/<code>GTR</code> are IQ-TREE's own fixed built-in equal-structure multistate
|
||
model, not a hook for an arbitrary custom matrix; a custom matrix is a
|
||
file path in <code>-m</code> directly, no <code>-mdef</code> or <code>GTRX</code> involved.)</em></p>
|
||
<p>Source for the empirical findings above: <code>model/modelmorphology.cpp</code> and
|
||
<code>model/modelmarkov.cpp</code> in the local <code>iqtree/iqtree3</code> source (cloned to
|
||
inspect the exact parsing logic after documentation didn't resolve the
|
||
<code>+Fname</code> reference error) — more reliable here than the PDF/web manual,
|
||
which (like TNT/PhyG) doesn't always match this specific binary. Original
|
||
(partially superseded) sources: <a href="https://iqtree.github.io/doc/Substitution-Models">Substitution
|
||
Models</a>, <a href="https://iqtree.github.io/doc/Complex-Models">Complex
|
||
Models</a>.</p>
|
||
<p><strong>Relation to the existing calibration.</strong> <code>sub_cost[a][b] = -ln(rate)</code>
|
||
(see "A concrete Sankoff cost matrix" above) is already a log-rate — a
|
||
genuine CTMC rate matrix <code>Q</code> could plausibly be recovered as
|
||
<code>rate(a,b) = exp(-cost(a,b))</code>, renormalised so each row sums to zero, once
|
||
the stationary-frequency gap above is closed. Superseded by the two
|
||
sections below, which implement and then substantially revise this.</p>
|
||
<h3 id="r-implemented-exclude-genome-added-rogue-taxon-test-negative-2026-08-12"><code>R</code>/<code>π</code> implemented; <code>--exclude-genome</code> added; rogue-taxon test negative (2026-08-12)</h3>
|
||
<p><code>write_iqtree</code> (<code>obikmer/src/cmd/phylo/iqtree.rs</code>) implements exactly
|
||
the <code>R = exp(-cost)</code> / empirical-<code>π</code> design above: writes
|
||
<code><prefix>_iqtree.model</code> (lower-triangular <code>R</code>, PAML order, then <code>π</code>) and
|
||
<code><prefix>_iqtree.fasta</code> (alignment recoded to the compact <code>0..k-1</code>
|
||
alphabet of states actually present), states/frequencies restricted to
|
||
whichever of the 16 canonical states actually occur, compactly
|
||
renumbered — see the <code>MORPH{N}</code> risk above. Verified end to end on the
|
||
20-genome benchmark: real, non-zero, varied branch lengths (<code>Total tree
|
||
length: 6.773</code>), converged log-likelihood, <code>State frequencies:</code> in
|
||
IQ-TREE's own output matching the computed <code>π</code> exactly.</p>
|
||
<p><strong>Rogue-taxon hypothesis, tested and refuted.</strong> The backbone (<em>Yersinia</em>,
|
||
<em>Proteus</em>, <em>Opitutus</em>, <em>Shouchella</em>, <em>Wolbachia</em>...) resolves as a
|
||
near-linear comb with ~8 near-zero branch lengths — visible on the real
|
||
tree, and independently on a Mash+NJ tree built from an entirely
|
||
different signal (whole-genome k-mer distance, no relation to the Sankoff
|
||
pipeline), which shows the same comb shape. Hypothesis: <code>Saccharolobus</code>/
|
||
<code>Candidozyma</code> (near-zero real signal — see below) destabilise the
|
||
heuristic tree search enough to also blur resolution elsewhere ("rogue
|
||
taxa", a documented phenomenon). Tested directly: reran <code>iqtree3</code> after
|
||
removing both taxa (<code>--exclude-genome</code>, <code>+ASC</code> recomputed on the
|
||
resulting variable-sites-only alignment since removing taxa turns some
|
||
columns invariant) — <strong>still exactly 8 near-zero backbone branches</strong>,
|
||
identical comb shape. Refuted for this dataset: the backbone's weak
|
||
resolution is a property of the character system's signal at that
|
||
divergence depth (matches "Run 3" above), not rogue-taxon interference.</p>
|
||
<p><strong>Diagnosed why <code>Saccharolobus</code>/<code>Candidozyma</code> place so poorly</strong>, using a
|
||
new diagnostic (<code>--raw-snp-counts</code>, <code><prefix>_rawsnp_counts.csv</code>: one row
|
||
per genome pair, <code>n_snp,n_shared,n_eligible,ratio</code> — the counts
|
||
<code>--raw-snp-distance</code>'s ratio-only matrix discards, needed because
|
||
<code>ratio=0.0</code> from 2 eligible loci and from 2000 look identical in the
|
||
ratio alone). <code>Saccharolobus</code> has 21,020 non-<code>∅</code> sites in the real
|
||
alignment (comparable to other taxa) — but <strong>21,017 of them (100.0%) are
|
||
private</strong>: no other genome has a non-<code>∅</code> state at the same site. Only 3
|
||
sites are shared with any other genome at all (1, 2, and 14
|
||
co-occurring genomes respectively). Real informativeness for placement
|
||
tracks shared sites, not raw non-<code>∅</code> count — with no other archaeon (or
|
||
even archaea-adjacent bacterium) in the dataset, there is nothing to
|
||
anchor <code>Saccharolobus</code>'s position to, regardless of how much of its own
|
||
data exists. <code>Candidozyma</code> shows the same pattern, more extreme (almost
|
||
all <code>NA</code>, the 3 non-<code>NA</code> pairs all exactly <code>0.0</code>, never <code>1.0</code> — itself a
|
||
tell: with <code>n_eligible=1</code>, the ratio can only be exactly <code>0</code> or <code>1</code>, so
|
||
3-for-3 landing on <code>0.0</code> is more than sampling noise alone would predict;
|
||
possibly ascertainment bias — see below — or possibly a few genuinely
|
||
ultra-conserved loci; not resolved).</p>
|
||
<p><strong><code>--exclude-genome LABEL</code></strong> (repeatable, <code>obikmer phylo</code>) added for
|
||
exactly this kind of test: zeroes the excluded genome's row/column in
|
||
<code>RawSnpDistanceOutput</code> after <code>raw_snp_distance</code> runs (a pair with zero
|
||
counts is already skipped by <code>calibrate_p_hat</code>/<code>base_pair_tally</code> — no
|
||
<code>obikindex</code> traversal change needed) and drops its row from
|
||
<code>SnpAlignment</code> before any output is written. Deliberately <em>not</em> index
|
||
surgery (a new, smaller on-disk index) — genome sets to exclude are
|
||
expected to change between quick tests, so an in-memory filter is the
|
||
right tool, not a new index-rewriting subsystem. Scoped to the sibling-annex
|
||
family of computations (<code>--raw-snp-distance</code>/<code>--raw-snp-counts</code>/<code>--snp</code>/
|
||
<code>--sankoff</code> and everything it implies) — does not affect the plain
|
||
<code>--metric</code> distance matrix/NJ/UPGMA path (a different, unrelated
|
||
computation on <code>idx.distance()</code>, not touched).</p>
|
||
<p><strong>Caveat surfaced while reusing a stale <code>π</code></strong>: rerunning IQ-TREE on a
|
||
genome-reduced alignment while keeping the <em>original</em> (20-genome)
|
||
model file is inconsistent — <code>π</code>'s composition shifts once low-cardinality-0
|
||
columns that were only variable because of the removed taxa drop out.
|
||
Measured directly: <code>π(∅)</code> <code>0.9048</code> (20 genomes) → <code>0.8869</code> (18 genomes,
|
||
variable sites only) — real (~378k affected cells, from <code>Saccharolobus</code>'s
|
||
~21k private sites × 18 remaining genomes) but modest in <em>proportion</em>
|
||
(~2 points) because the alignment was already monomorphic-filtered
|
||
before <code>Saccharolobus</code> was ever added, so removing it only drops the
|
||
subset of columns that were variable <em>because of</em> it specifically, not
|
||
every column it appears in.</p>
|
||
<h3 id="r-via-exp-cost-is-wrong-for-a-ctmc-cardinalitycomposition-decomposition-open-2026-08-12"><code>R</code> via <code>exp(-cost)</code> is wrong for a CTMC; cardinality/composition decomposition (open, 2026-08-12)</h3>
|
||
<p><strong>The flaw in <code>R = exp(-cost)</code>, precisely.</strong> <code>cost</code> (<code>build_cost_matrix</code>'s
|
||
output) is a <em>shortest-path closure</em> over an elementary-edit graph
|
||
(Floyd-Warshall) — correct and required for Sankoff parsimony, where
|
||
<code>cost(a,b)</code> must be a metric. But a CTMC's own matrix exponential
|
||
(<code>exp(Qt)</code>, computed internally by IQ-TREE) <em>already</em> sums over
|
||
paths of every length through the elementary rates — that's the whole
|
||
mechanism by which a CTMC generates indirect transitions. Feeding it a
|
||
pre-summed, multi-hop shortest-path cost and exponentiating that <em>again</em>
|
||
as if each entry were a direct edge double-applies the "compose multiple
|
||
steps" logic once in log-space (Floyd-Warshall, additive) and once more
|
||
inside IQ-TREE's own exponential — systematically over-penalising
|
||
non-adjacent state pairs (e.g. <code>∅→{A,C}</code> priced as two chained edges,
|
||
<code>2×c_ctx</code>, when it should be one direct lookup).</p>
|
||
<p><strong>Resolution, in two parts — both estimated directly from the real
|
||
alignment, not smoothed through a small parametric formula:</strong></p>
|
||
<ol>
|
||
<li><strong>Cardinality model</strong>: a 5-state (<code>0,1,2,3,4</code>) first-order Markov
|
||
chain, estimated from the empirical cardinality co-occurrence table
|
||
(pooled across all <em>included</em> genome pairs — same
|
||
saturated/<code>NA</code>-pair exclusion discipline as <code>calibrate_p_hat</code>, 14
|
||
saturated + 18 <code>NA</code> pairs excluded of 190 in the benchmark; the
|
||
result is materially different from the unfiltered version and more
|
||
internally consistent, not just "cleaner"). Diagonal included (the
|
||
probability of a family <em>staying</em> at the same cardinality is part of
|
||
the model, not assumed away). Measured, cost <code>= -ln(observed/expected
|
||
under independence)</code>, on the 20-genome benchmark (<code>π</code> here from the
|
||
marginal cardinality distribution: <code>90.481%, 5.261%, 4.183%, 0.067%,
|
||
0.008%</code> for <code>c=0..4</code> respectively — note <code>c=2</code> is <em>not</em> rare, almost
|
||
as common as <code>c=1</code>):</li>
|
||
</ol>
|
||
<p><div class="highlight"><pre><span></span><code>c=0/c=0: 0.166 c=0/c=1: 0.160 c=0/c=2: 0.157 c=0/c=3: 0.059 c=0/c=4: 0.006
|
||
c=1/c=1: -0.803 (enriched) c=1/c=2: 0.381 c=1/c=3: 1.027 c=1/c=4: 1.794
|
||
c=2/c=2: 4.032 (sharply suppressed) c=2/c=3: 2.620 c=2/c=4: 2.416
|
||
c=3/c=3: -0.254 c=3/c=4: -0.411 (too few observations to trust)
|
||
</code></pre></div>
|
||
Two robust findings, both stable under the saturation filter: (a) all
|
||
<code>∅</code>-involving costs are low and close to each other (<code>0.006</code>–<code>0.166</code>)
|
||
regardless of how many members are gained/lost at once — sharply at
|
||
odds with the current model's implicit <code>×2</code>/<code>×3</code> multi-hop scaling;
|
||
(b) <code>c=2/c=2</code> (two genomes both showing an ambiguous 2-member state at
|
||
the same site) is dramatically under-represented (~1.5–1.8% of the
|
||
independence expectation) — a real, robust anomaly, not explained.</p>
|
||
<ol>
|
||
<li><strong>Composition model</strong>: unchanged — <code>sub_cost</code> (the existing 6-category
|
||
Ts/Tv-biased calibration), estimated <strong>only from unambiguous sites</strong>
|
||
(cardinality-1 ↔ cardinality-1 pairs), because composition bias is a
|
||
substitution phenomenon and only means something when cardinality is
|
||
conserved. Checked whether composition bias also appears in pure
|
||
gain/loss events (no substitution involved, so no bias expected a
|
||
priori): single-base "gain" (<code>∅→</code>single base) is close to uniform
|
||
(<code>24.35/25.92/26.22/23.51%</code> for A/C/G/T) — consistent with "no
|
||
mutational mechanism, so no bias" as expected. Two-base "gain"
|
||
(<code>∅→</code>2-member state) is <em>not</em> uniform even after correcting for the
|
||
real (non-25/25/25/25) marginal base frequencies: <code>{A,G}</code>/<code>{C,T}</code>
|
||
(the transition-linked pairs) mildly enriched (<code>obs/exp</code> <code>1.12</code>/
|
||
<code>1.16</code>), <code>{C,G}</code> sharply suppressed (<code>obs/exp 0.563</code>, expected to be
|
||
the <em>most</em> common pair under independence since C and G are
|
||
individually the two most frequent bases, observed the least) — an
|
||
unexplained anomaly, deliberately <strong>not</strong> built into the model (no
|
||
mechanistic story for why gain/loss would carry a <code>{C,G}</code>-specific
|
||
bias), left as an open puzzle rather than fit.</li>
|
||
</ol>
|
||
<p><strong>Composing the two into a full pairwise cost — the actual replacement
|
||
for both the elementary-edge graph and its Floyd-Warshall closure.</strong> For
|
||
any two states <code>A</code>, <code>B</code> (not just the "clean" same-cardinality or
|
||
pure-subset cases the current graph handles directly): let <code>shared = A∩B</code>
|
||
(free), <code>lost = A\B</code>, <code>gained = B\A</code>. Pair off <code>k = min(|lost|,|gained|)</code>
|
||
elements between <code>lost</code> and <code>gained</code> as substitution events, choosing the
|
||
pairing that minimises total <code>sub_cost</code> (a trivial assignment problem —
|
||
at most 4 elements per side). The <code>|lost|-k</code> (or <code>|gained|-k</code>) leftover,
|
||
unpaired elements are a <em>pure</em> cardinality change, <code>|A|→|B|</code>, costed by
|
||
<strong>one direct lookup</strong> in the cardinality model above — no chaining, no
|
||
Floyd-Warshall. Example: <code>{A,C}→{G}</code> (cardinality 2→1, no shared base):
|
||
pair 1 substitution (cheaper of <code>A→G</code>, <code>C→G</code>), 1 base left over unpaired
|
||
→ cost <code>= sub_cost(chosen pair) + cardinality_cost(2→1)</code>. This is a
|
||
direct, closed-form cost for <em>every</em> pair of the 16 states, replacing
|
||
<code>build_cost_matrix</code>'s graph-plus-shortest-path construction outright —
|
||
and specifically fixes the CTMC double-counting problem, since every
|
||
entry is now a single decomposed lookup, never a sum of chained edges.</p>
|
||
<p><strong>Reframed as a likelihood (product of probabilities), not a cost (sum of
|
||
<code>-ln</code>s) — same content, but forces the diagonals (no-change cases) to be
|
||
kept rather than implicitly dropped.</strong> Both sub-models are proper
|
||
transition <em>probability</em> matrices, diagonal included — <code>P_cardinality</code>
|
||
includes "stay at the same cardinality", <code>P_composition</code> includes "stay
|
||
the same base" (e.g. <code>P(G→G)</code>, not just <code>sub_cost</code>'s off-diagonal
|
||
entries). Composing:</p>
|
||
<div class="highlight"><pre><span></span><code>P(A→B) = P_cardinality(|A|→|B|)
|
||
× ∏_{x ∈ A∩B} P_composition(x→x) (shared bases: "stayed")
|
||
× ∏_{(x,y) paired} P_composition(x→y) (parsimony-paired substitutions)
|
||
</code></pre></div>
|
||
<p>with the unpaired leftover <code>lost</code>/<code>gained</code> elements (if <code>|lost|≠|gained|</code>)
|
||
contributing <em>nothing further</em> beyond the <code>P_cardinality</code> term already
|
||
counted — consistent with the finding above that pure gain/loss carries
|
||
no separate composition bias worth modelling. Worked example,
|
||
<code>{A,C}→{A,G}</code> (<code>shared={A}</code>, one paired substitution <code>C→G</code>, nothing left
|
||
over): <code>P = P_cardinality(2→2) · P_composition(A→A) · P_composition(C→G)</code>.</p>
|
||
<p>The row-wise product of these two independently-calibrated models isn't
|
||
guaranteed to already sum to exactly <code>1</code> across all <code>B</code> for a fixed <code>A</code>
|
||
(the two aren't perfectly independent in reality) — so each row of the
|
||
resulting 16×16 matrix is renormalised (divided by its own sum) after
|
||
composition, <strong>not</strong> the matrix as a whole (which would produce a joint
|
||
distribution over <code>(A,B)</code> pairs, the wrong object — a transition matrix
|
||
needs each row, "given I start in <code>A</code>", to be a valid distribution over
|
||
where I end up).</p>
|
||
<p>This now covers every pair of the 16 states with no unhandled case
|
||
identified. Gives 120 parameters, but derived from two small,
|
||
well-estimated pieces (a 5×5 cardinality model, a 4×4 composition model)
|
||
rather than fit or smoothed independently per pair.</p>
|
||
<h3 id="implemented-2026-08-12-pairwise_cost_matrix-replaces-build_cost_matrix-entirely">Implemented (2026-08-12): <code>pairwise_cost_matrix</code> replaces <code>build_cost_matrix</code> entirely</h3>
|
||
<p>New module <code>obikindex/src/cardcomp.rs</code>, replacing <code>sankoff::build_cost_matrix</code>
|
||
and the <code>c_ctx</code>/<code>SankoffWeights</code>/<code>PHatEstimate</code>/<code>calibrate_p_hat</code>/
|
||
<code>c_ctx_from_p_hat</code>/<code>mean_substitution_cost</code>/<code>substitution_costs_from_tally</code>
|
||
machinery it depended on outright — not kept in parallel as a fallback
|
||
(all now unreferenced outside their own tests; <code>sankoff.rs</code> itself is a
|
||
pending removal, not yet done).</p>
|
||
<p><strong>New primitives, <code>obikindex/src/siblings.rs</code>:</strong>
|
||
- <code>BasePairTally</code> gained a <code>same: [u64; 4]</code> field (diagonal — "both
|
||
genomes at the same single base", pooled from
|
||
<a href="obikindex::KmerIndex::base_pair_tally"><code>base_pair_tally</code></a>'s existing
|
||
traversal, extended to also tally the <code>bi == bj</code> case it previously
|
||
discarded).
|
||
- <code>CardinalityTally { counts: [[u64; 5]; 5] }</code> and
|
||
<code>KmerIndex::cardinality_tally</code>, a new traversal (same shape as
|
||
<code>snp_pseudo_alignment</code>'s — needs full per-genome presence masks, not
|
||
<code>scan_family_pairs</code>'s single-resolved-form view, since cardinality 2-4
|
||
is exactly the signal being tallied, not noise to drop). Same
|
||
saturated/no-data pair exclusion as <code>base_pair_tally</code>. Restricted to
|
||
variable families (<code>family_size() >= 2</code>), matching
|
||
<code>snp_pseudo_alignment</code>'s own scope. Verified against the 20-genome
|
||
benchmark: counts match the earlier hand-rolled Python analysis exactly
|
||
(e.g. <code>c=0/c=0: 117,158,166</code>, <code>c=0/c=1: 13,707,223</code> — the same numbers
|
||
this whole investigation started from).</p>
|
||
<p><strong><code>cardcomp.rs</code>:</strong>
|
||
- <code>cardinality_transition_probs</code>/<code>composition_transition_probs</code>: row-
|
||
normalise the two tallies into proper transition probability matrices,
|
||
diagonal included ("stay the same" is a real, calibrated outcome).
|
||
- <code>pairwise_cost_matrix</code>: for every pair of the 16 states, <code>shared = A∩B</code>
|
||
contributes <code>∏ P_composition(x→x)</code>; <code>lost = A\B</code>, <code>gained = B\A</code> are
|
||
parsimony-paired (<code>best_pairing_cost</code>, brute-force over the ≤4!
|
||
injections — small enough that hand-rolling beats a dependency) into
|
||
substitution events on <code>P_composition</code>, minimising total <code>-ln</code>; the
|
||
cardinality-difference leftover is priced once via
|
||
<code>P_cardinality(|A|→|B|)</code>, never chained. Row-normalised, <code>-ln</code>'d, then
|
||
<strong>symmetrised</strong>: <code>cost_sym(A,B) = (cost(A,B)+cost(B,A))/2</code> — equivalent
|
||
to taking the <em>geometric</em> mean of the two raw probabilities
|
||
(<code>-ln(√(P(A,B)·P(B,A))) = (-ln P(A,B) - ln P(B,A))/2</code>), not their
|
||
arithmetic mean. Required, not just convenient for IQ-TREE's
|
||
lower-triangular file format: Sankoff parsimony's score is independent
|
||
of where an <em>unrooted</em> tree (what TNT/PhyG actually search over) gets
|
||
rooted only if the cost matrix is symmetric — the discrete-parsimony
|
||
analogue of CTMC reversibility, established by direct reasoning, not
|
||
assumed. Bonus of the same decision: 120 free parameters instead of the
|
||
240 a fully asymmetric matrix would need.</p>
|
||
<p><strong>Verified on the 20-genome benchmark</strong>: resulting matrix symmetric
|
||
(checked numerically, zero asymmetric cells), zero diagonal, no NaN/Inf.
|
||
<code>--tnt</code> output still loads into TNT with no triangle-inequality warning
|
||
(<code>scaled_metric_matrix</code>'s rounding-metric-closure step still needed and
|
||
still applied — nothing in the new construction guarantees the <em>rounded
|
||
integer</em> matrix stays a metric, even though the real-valued one is exact
|
||
by construction here, unlike the old Floyd-Warshall-closed matrix which
|
||
needed it for a different reason). IQ-TREE loads the new model file and
|
||
reports the same <code>π</code> as before (only <code>R</code> changed).</p>
|
||
<h3 id="two-consistency-bugs-found-and-fixed-post-implementation-2026-08-13">Two consistency bugs found and fixed post-implementation (2026-08-13)</h3>
|
||
<p><strong><code>--exclude-genome</code> didn't drop columns that become monomorphic once the
|
||
excluded genome(s) are gone.</strong> <code>snp_pseudo_alignment</code>'s "variable family"
|
||
test (<code>family_size() >= 2</code>) is a property of the annex computed over
|
||
<em>every</em> genome in the index — unaffected by the CLI-level exclusion, which
|
||
only dropped the excluded genome's <em>row</em>. A family variable only because
|
||
of the excluded genome stayed in the alignment as a now-constant column —
|
||
silently wrong data for TNT/PhyG, and a hard failure for IQ-TREE's <code>+ASC</code>
|
||
(verified: excluding 2 taxa on the benchmark left 116,351 such columns —
|
||
matches the manual <code>+ASC</code> failures hit earlier in this same investigation,
|
||
before <code>--exclude-genome</code> existed). Fixed in <code>drop_excluded</code>
|
||
(<code>obikmer/src/cmd/phylo/mod.rs</code>): after dropping excluded rows,
|
||
re-scan each column among the <em>surviving</em> sequences and drop any that are
|
||
now constant. Verified: 908,723 → 792,372 sites after excluding 2 taxa,
|
||
zero monomorphic columns remain, <code>π</code> recomputed from the corrected
|
||
alignment matches an independent recount exactly. The compact-alphabet
|
||
renumbering (<code>iqtree::compact_alphabet</code>) needed no equivalent fix — it
|
||
already recomputes which of the 16 states occur fresh on every call, from
|
||
whatever alignment it's actually handed, so a symbol disappearing (e.g.
|
||
excluding every genome that carries <code>N</code>) is already handled correctly;
|
||
verified directly (excluded 13 genomes to force <code>N</code> out: "15 of 16 states"
|
||
reported, correctly-shaped model file).</p>
|
||
<p><strong><code>cardinality_tally</code>'s <code>family_size() >= 2</code> filter looked inconsistent
|
||
with <code>base_pair_tally</code> — removing it was tried, and was wrong; reverted.</strong>
|
||
<code>cardinality_tally</code> (modelled after <code>snp_pseudo_alignment</code>) had the
|
||
filter; <code>base_pair_tally</code> didn't (it visits every family via
|
||
<code>scan_family_pairs</code> unconditionally, folding fully-invariant loci into its
|
||
own <code>same</code> diagonal). Read as <code>cardinality_tally</code> under-counting its
|
||
diagonal relative to <code>base_pair_tally</code>, and — independently — as another
|
||
angle on the <code>--exclude-genome</code> drift (<code>family_size()</code> being global-only
|
||
meant a family kept here post-exclusion could differ from what the
|
||
now-correctly-filtered alignment kept). First fix tried: drop
|
||
<code>cardinality_tally</code>'s filter entirely, matching <code>base_pair_tally</code>'s
|
||
whole-annex scope.</p>
|
||
<p><strong>That fix was empirically wrong, confirmed by a real IQ-TREE run, not
|
||
just a hunch.</strong> Log-likelihood dropped from the earlier correct run's
|
||
<code>-8,364,671</code>/<code>-8,371,082</code> to <code>-9,170,228</code> (worse fit, not better), with
|
||
repeated <code>NNI search needs unusual large number of steps (20) to
|
||
converge!</code> warnings — and the completed run's <strong>total tree length came
|
||
out at 67.644</strong>, roughly 30× the earlier correct runs' ~2.0, i.e. branches
|
||
blowing up/saturating. Root cause, only clear in hindsight: <code>+ASC</code>
|
||
("ascertainment bias correction") exists specifically because the
|
||
likelihood only ever sees <em>variable</em> sites — the alignment fed to
|
||
IQ-TREE, by construction, contains not one invariant column. Calibrating
|
||
<code>R</code> from a population overwhelmingly dominated by genome-wide invariant
|
||
background (family_size()<2 loci outnumber the ~908k variable ones by
|
||
orders of magnitude) describes a completely different population than the
|
||
one <code>+ASC</code> and the alignment actually model — the "consistency" argument
|
||
for matching <code>base_pair_tally</code>'s scope was real, but pointed the wrong
|
||
way: <code>base_pair_tally</code>'s own unrestricted <code>same</code> diagonal turned out to
|
||
have the <em>identical</em> latent bug (only unmasked once its diagonal existed
|
||
at all, which happened earlier the same day when <code>same</code> was added), not a
|
||
correct baseline to match <code>cardinality_tally</code> to.</p>
|
||
<p><strong>Final fix</strong>: restored <code>cardinality_tally</code>'s <code>family_size() >= 2</code> filter,
|
||
and gave <code>base_pair_tally</code>'s <code>same</code> diagonal the equivalent restriction —
|
||
<code>scan_family_pairs</code> (shared with <code>raw_snp_distance</code>, which legitimately
|
||
<em>does</em> want fully-invariant families counted as <code>shared</code>) now passes an
|
||
extra <code>variable: bool</code> (the family's own <code>family_size() >= 2</code>) to its
|
||
<code>on_pair</code> callback; <code>base_pair_tally</code> only increments <code>same</code> when
|
||
<code>variable</code> is true, <code>raw_snp_distance</code>'s callback ignores the new
|
||
argument. Both tallies now describe the same variable-families-only
|
||
population the <code>+ASC</code>-corrected alignment does. Verified: calibration
|
||
counts back to their original values exactly (<code>c=0/c=0</code>: 117,158,166,
|
||
matching the pre-regression run bit for bit), and a full IQ-TREE rerun
|
||
converged normally — log-likelihood <code>-8,389,106.273</code> (same order as the
|
||
two earlier correct runs), <strong>total tree length 2.044</strong> (was 67.644), no
|
||
NNI convergence warnings.</p>
|
||
<h2 id="heterozygosity-ploidy-and-consensus-assembly-inputs">Heterozygosity, ploidy, and consensus-assembly inputs</h2>
|
||
<p>A within-genome multiplicity signal (more than one of the 4 central forms
|
||
present at a locus) is produced identically by two distinct causes:
|
||
paralogous duplication and diploid/polyploid heterozygosity. K-mer data alone
|
||
cannot distinguish them. The one real discriminator is sequencing depth
|
||
(heterozygous site: total depth of the present forms ~= the genome's
|
||
single-copy average; duplication: ~2x or more) — but that signal only exists
|
||
if genome "counts" are raw-read depth (FASTQ input), not occurrence counts in
|
||
an assembled FASTA, where per-locus depth is not preserved.</p>
|
||
<p><strong>Magnitude is taxon- and mating-system-dependent, not universal.</strong>
|
||
Heterozygosity density: mammals ~1 site / 1-1.5 kb (~0.1%); highly
|
||
outcrossing plants (maize, poplar) reported an order of magnitude higher
|
||
(~1%); self-fertilising plants (<em>Arabidopsis thaliana</em>) near zero — but with
|
||
a documented failure mode where segmental duplication masquerades as
|
||
"pseudo-heterozygosity"; fungi split between haploid vegetative stages
|
||
(non-issue) and dikaryotic Basidiomycetes, where two long-diverged haploid
|
||
nuclei coexist without fusing. The estimator's target use case (closely
|
||
related genomes, k=31) is exactly where the stringent filter above costs the
|
||
least for low-heterozygosity taxa and the most for outcrossing/dikaryotic
|
||
ones — no universal threshold; this is a scope caveat to document, not a
|
||
problem to solve generically.</p>
|
||
<p><strong>Why assembled-consensus inputs don't make measured distances wrong.</strong>
|
||
Phylogenetic inputs are near-universally assemblies, not raw reads, and
|
||
assemblers collapse heterozygous sites to one consensus allele per
|
||
position — effectively an arbitrary, largely uncorrelated-between-assemblies
|
||
choice at each het site. This does not inject unbounded noise: standard
|
||
population genetics gives <code>d_xy = d_a + (pi_A + pi_B)/2</code> — the expected
|
||
pairwise difference between a random allele of population A and a random
|
||
allele of population B equals the net (fixed) divergence <code>d_a</code> plus the
|
||
average of the two populations' own within-population diversity <code>pi</code>.
|
||
Consensus flattening realises exactly this random-allele draw, so the
|
||
measured genome-to-genome distance is a <code>d_xy</code>-like quantity, not <code>d_a</code> —
|
||
inflated by heterozygosity by a well-characterised additive term, not
|
||
distorted unpredictably. The term is negligible when <code>pi << d_xy</code> (the common
|
||
case for cross-species comparisons), and becomes material precisely in the
|
||
two cases already flagged above: very closely related genomes (this
|
||
estimator's explicit target) and highly heterozygous outcrossing organisms,
|
||
where <code>pi</code> and <code>d_xy</code> are the same order of magnitude.</p>
|
||
<p>Caveat: this assumes the flattening is uncorrelated with the phylogenetic
|
||
signal — plausible for de novo assembly, not guaranteed for reference-guided
|
||
assembly biased toward one allele (e.g. the reference's) at each het site,
|
||
which would turn the noise term into a systematic bias toward the reference
|
||
lineage. Not evaluated here.</p>
|
||
<p><strong>Forward-looking implication, not part of the current design.</strong> The
|
||
multiplicity > 1 signal discarded by the stringent filter is a crude
|
||
per-genome proxy for <code>pi</code> (under low background paralogy). If a <code>pi_hat</code> per
|
||
genome were tallied alongside <code>SnpTally</code>, a <code>d_a</code> correction
|
||
(<code>p_hat - mean(pi_hat_i, pi_hat_j)/2</code>, roughly) could recover an estimate
|
||
closer to net divergence instead of <code>d_xy</code> — a possible extension, not
|
||
scoped here.</p>
|
||
<h2 id="sufficient-statistic-4x4-base-pair-tally">Sufficient statistic: 4x4 base-pair tally</h2>
|
||
<p>Tabulating the joint distribution of <code>(center_i, center_j)</code> over conserved-flank
|
||
loci, per genome pair, is sufficient for every downstream correction:</p>
|
||
<table>
|
||
<thead>
|
||
<tr>
|
||
<th>Estimator</th>
|
||
<th>Input</th>
|
||
<th>Formula</th>
|
||
</tr>
|
||
</thead>
|
||
<tbody>
|
||
<tr>
|
||
<td>Raw p-distance</td>
|
||
<td>total off-diagonal / total</td>
|
||
<td><code>p = SNP / (SNP + shared)</code></td>
|
||
</tr>
|
||
<tr>
|
||
<td>Jukes-Cantor</td>
|
||
<td>p</td>
|
||
<td><code>d = -3/4 * ln(1 - 4p/3)</code></td>
|
||
</tr>
|
||
<tr>
|
||
<td>Kimura 2-parameter</td>
|
||
<td>transition rate P, transversion rate Q</td>
|
||
<td><code>d = 1/2 ln(1/(1-2P-Q)) + 1/4 ln(1/(1-2Q))</code></td>
|
||
</tr>
|
||
<tr>
|
||
<td>LogDet/paralinear</td>
|
||
<td>full 4x4 + base-composition margins</td>
|
||
<td><code>d ~= -1/4 ln det(F)</code>, robust to non-stationary base composition</td>
|
||
</tr>
|
||
</tbody>
|
||
</table>
|
||
<p>JC/K2P need only the total and the transition/transversion split (the
|
||
diagonal collapses to a single "shared" total). LogDet needs the full 4x4,
|
||
already populated at no extra cost (see Step 1/2 below).</p>
|
||
<p>Memory for the 4x4 tally: <code>n^2 * 16</code> counters. Trivial for the project's
|
||
genome-scale use case (tens to hundreds of genomes); ~13 GB at n=10^4 — outside
|
||
scope but worth flagging if n grows.</p>
|
||
<h2 id="biases-properties-of-the-estimator-not-defects">Biases (properties of the estimator, not defects)</h2>
|
||
<ol>
|
||
<li><strong>Conserved-flank ascertainment bias.</strong> Only SNPs with intact <code>2m</code>-base
|
||
flanks are visible; window-intact probability decays as <code>(1-p)^{2m}</code>. For
|
||
k=31 (2m=30): 0.74 at p=1%, 0.21 at p=5%, 0.04 at p=10%. This estimator
|
||
targets <strong>closely related genomes</strong>. Under rate heterogeneity across sites
|
||
(universal in practice), conserved flanks correlate with slow centers, so
|
||
<code>p_hat</code> underestimates the genome-wide average rate — it specifically
|
||
estimates the substitution rate of <strong>conserved regions</strong>.
|
||
Two distinct factors are at play here, not one: <code>P(centre of a given
|
||
window is a SNP) = p</code> exactly, <strong>independent of k</strong> — a direct restatement
|
||
of the raw per-site rate via the bijective window<->centre-position
|
||
correspondence (Statistic section above), not a k-dependent quantity.
|
||
<code>(1-p)^{2m}</code> is the <em>separate</em>, genuinely k-dependent ascertainment factor
|
||
(are the flanks also intact). The two multiply:
|
||
<code>P(usable window showing a central SNP) = p * (1-p)^{2m}</code> — e.g. at
|
||
p=1/31 (~3.2%), k=31: <code>p * (1-p)^30 ~= 0.0323 * 0.374 ~= 1.2%</code>, i.e. about
|
||
1 window in 83, not 1 in 31 (which is only the centre-mutated fraction,
|
||
before requiring intact flanks).</li>
|
||
<li><strong>Bias toward isolated SNPs.</strong> Two SNPs within k of each other disqualify
|
||
each other's flanks. Hypervariable regions are invisible by construction.</li>
|
||
<li><strong>Indels are invisible.</strong> A frameshift destroys k-mer matches in a block;
|
||
this channel captures substitutions only. Indel divergence shows up as lost
|
||
shared k-mers (lower Jaccard/Mash), not as SNP signal.</li>
|
||
<li><strong>k-dependent specificity.</strong> "A k-mer match implies common ancestry" is
|
||
quantitative. For a 3 Gbp genome, expected random flank-30 collisions
|
||
(k=31): <code>(3e9)^2 / 4^30 ~= 8</code> — negligible. At k=21: <code>(3e9)^2 / 4^20 ~= 2e6</code>
|
||
— no longer negligible. k=31 is safe; k<=21 is marginal to unreliable for
|
||
large genomes. The large k that guarantees homology is the same k that
|
||
shrinks the detectable-divergence window — an inherent tension.</li>
|
||
</ol>
|
||
<h2 id="implementation-avoid-materializing-a-de-bruijn-graph">Implementation: avoid materializing a de Bruijn graph</h2>
|
||
<p>A central-SNP pair is topologically a simple bubble in the colored de Bruijn
|
||
graph (source/sink k-mer shared, two length-k branches differing only at the
|
||
midpoint). Classical bubble-calling (Cortex/discoSNP-style) finds these, but
|
||
requires the graph — nodes plus adjacency for ~10^9 colored k-mers — resident
|
||
in memory. <strong>Rejected</strong>: prohibitive RAM for this project's scale.</p>
|
||
<p>A naive per-pair generalisation of variant lookup across n genomes (query each
|
||
non-shared k-mer's 3 central variants against every counterpart genome's
|
||
index) costs <code>O(n^2 . N . 3)</code> random lookups, with the same k-mer's 3 variants
|
||
regenerated and requeried once per counterpart genome — pure redundant work.
|
||
<strong>Rejected</strong> as the basis for an n-genome design.</p>
|
||
<h2 id="implementation-sequential-per-partition-sweep-no-scratch-no-graph">Implementation: sequential per-partition sweep (no scratch, no graph)</h2>
|
||
<p><code>KmerIndex::distance()</code> already opens every partition's <code>presence_store</code>/
|
||
<code>count_store</code> simultaneously, memory-mapped, into one <code>LayeredStore</code>
|
||
(<code>distance.rs:73-77</code>). "Querying another partition" is therefore not a new
|
||
I/O pattern to design — it is the same O(1) MPHF+evidence lookup the <code>query</code>
|
||
command already performs at scale. This lets the SNP tally be computed with
|
||
<strong>no scratch files and no auxiliary graph</strong>, by sweeping partitions once each
|
||
as a source:</p>
|
||
<ol>
|
||
<li>For source partition <code>p</code>, enumerate its <strong>distinct</strong> k-mers (one per MPHF
|
||
slot; each already carries its full multi-genome presence/count vector —
|
||
no need to explode per (k-mer, genome) occurrence).</li>
|
||
<li>For each, generate the 3 central-substitution variants and <strong>canonicalise
|
||
each independently</strong> (<code>min(kmer, revcomp)</code>, exactly as any normal query) —
|
||
this avoids the orientation edge case a masked-flank grouping would have
|
||
(a substitution that flips canonical orientation is handled correctly
|
||
because each variant is canonicalised on its own, not inferred from a
|
||
fixed-orientation flank key).</li>
|
||
<li>Compute each variant's target partition <code>q</code> via its minimizer; batch/sort
|
||
the partition's outgoing variant queries by <code>q</code> for locality.</li>
|
||
<li>Look up each variant in <code>q</code>'s already-mmap'd MPHF+evidence; on a hit,
|
||
combine the source's presence vector (base <code>a</code>) with the variant's
|
||
presence vector (base <code>b</code>): for every <code>i</code> carrying <code>a</code> and <code>j</code> carrying
|
||
<code>b</code>, <code>tally[i,j][a,b] += 1</code>.</li>
|
||
</ol>
|
||
<p><strong>Deduplication needs no persisted state.</strong> Sweeping partitions in a fixed
|
||
order <code>p = 0, 1, ..., P-1</code> and only acting on a variant when its target
|
||
partition <code>q >= p</code> guarantees each unordered SNP pair is counted exactly
|
||
once: a pair with <code>q < p</code> was already resolved earlier, when <code>q</code> was itself
|
||
the source partition and <code>p</code> (being <code>>= q</code>) was a valid forward target. No
|
||
cross-partition flag array is needed — the sweep order <em>is</em> the
|
||
deduplication rule. Within the same partition (<code>q == p</code>), a lightweight
|
||
transient tie-break suffices: either a <code>#slots(p)</code>-bit scratch flag reset per
|
||
partition, or simply comparing the two k-mers' raw <code>u64</code> encodings and only
|
||
counting when <code>kmer_source < kmer_variant</code> — no storage at all.</p>
|
||
<p>This is a distinct computation stage, not a <code>partial_*</code> in the existing
|
||
additive-by-partition sense: step 3-4 read across partition boundaries by
|
||
construction, unlike the row-local <code>partial_jaccard</code>/<code>partial_threshold_jaccard</code>
|
||
primitives. But it requires no new index files, permanent or scratch:
|
||
<code>unitigs.bin</code>, <code>mphf.bin</code>, <code>evidence.bin</code>, and the presence/count columns are
|
||
read as-is, and the only extra memory is the current partition's small
|
||
outgoing-query batch (<code>#kmers(p) * 3</code>, released once <code>p</code> is done) plus the
|
||
persistent <code>tally</code> accumulator (<code>n^2 * 16</code> counters, see above).</p>
|
||
<p><strong>Outer loop (over source partitions <code>p</code>) must stay sequential.</strong> Two
|
||
independent reasons, not just one: (a) memory — the bounded-footprint claim
|
||
above only holds with one partition's outgoing-query batch in flight; running
|
||
<code>T</code> source partitions concurrently multiplies that batch by <code>T</code>, exactly the
|
||
blowup the design avoids; (b) correctness — the <code>q >= p</code> deduplication rule
|
||
requires partitions to be claimed as sources in a fixed order; running <code>p1 <
|
||
p2</code> concurrently gives no guarantee <code>p1</code> has finished claiming its <code>q >= p1</code>
|
||
targets before <code>p2</code> starts claiming its own, breaking the "counted exactly
|
||
once" property.</p>
|
||
<p><strong>Inner loop (target-partition lookups for a fixed <code>p</code>) parallelises safely.</strong>
|
||
Each lookup is an O(1) read against an already-mmap'd structure, independent
|
||
of the others, with no growing allocation — no memory blowup, no ordering
|
||
dependency between different <code>q</code>. The only shared mutable state is <code>tally</code>;
|
||
give each worker thread a <strong>thread-local partial tally</strong> (fixed <code>n^2 * 16</code>
|
||
size, independent of partition size) and merge into the global <code>tally</code> once
|
||
<code>p</code>'s inner loop completes — the same reduce-then-merge pattern Rayon already
|
||
uses elsewhere in this codebase to open partitions in parallel. Extra memory:
|
||
<code>#threads * n^2 * 16</code>, negligible (~512 MB at n~1000, 32 threads) and
|
||
unrelated to partition size.</p>
|
||
<p><strong>Cost</strong>: <code>3 * N_distinct</code> MPHF lookups total across the whole index (each
|
||
partition swept once as source) — the same order of magnitude and the same
|
||
operation as running <code>query</code> over the index's entire k-mer content against
|
||
itself, three times. This is the tool's already-optimized regime, not a new
|
||
I/O profile to validate.</p>
|
||
<h2 id="cheaper-subsampling">Cheaper: subsampling</h2>
|
||
<p>Since the target is a ratio, restricting the source-partition sweep to a
|
||
bottom-<code>s</code> hash sketch (only enumerate k-mers with <code>hash < threshold</code> as
|
||
sources) divides the lookup count by the sampling factor without biasing
|
||
<code>p_hat</code>. Mash-like tradeoff: rate estimated from a sample, not the full
|
||
k-mer set.</p>
|
||
<h2 id="recommendation">Recommendation</h2>
|
||
<p>Sequential per-partition sweep (Route D): reuse the already-mmap'd
|
||
per-partition MPHF/evidence/presence structures for O(1) variant lookups,
|
||
dedup via the fixed sweep-order rule (<code>q >= p</code>, plus an in-partition
|
||
tie-break), no scratch files, no graph materialisation. Both the SNP
|
||
(off-diagonal) and shared (diagonal) counts are accumulated by this same
|
||
sweep, under the locus-eligibility rule chosen (raw or paralogy-filtered) —
|
||
not reused from the general-purpose <code>shared_kmers</code> matrix, whose raw-identity
|
||
definition does not apply the same copy-number constraint. Distances (p, JC,
|
||
K2P, LogDet) as finalisations of the resulting 4x4 tally, mirroring the
|
||
<code>partial_* -> *_dist_matrix</code> pattern used for Jaccard/Mash/Bray-Curtis/etc.</p>
|
||
<h2 id="detailed-implementation-plan">Detailed implementation plan</h2>
|
||
<p>Grounded in the current codebase. File/type references are anchors, not
|
||
prescriptions; adjust to reality when implementing.</p>
|
||
<h3 id="step-0-new-low-level-primitives-obikseq">Step 0 — new low-level primitives (<code>obikseq</code>)</h3>
|
||
<p>Two helpers do not yet exist and are prerequisites:</p>
|
||
<ol>
|
||
<li><strong>Central neighbours.</strong> <code>CanonicalKmerOf<L></code> already exposes
|
||
<code>left_canonical_neighbors()</code> / <code>right_canonical_neighbors()</code>
|
||
(<code>obikseq/src/kmer.rs</code>), each returning the 4 canonicalised neighbours at
|
||
an end position. Add <code>central_canonical_neighbors()</code> returning the 4
|
||
variants at position <code>m = (k-1)/2</code> (each independently canonicalised via
|
||
<code>.canonical()</code>). The 3 that differ from the source are the query variants;
|
||
skip the identity. Building on <code>nucleotide(i)</code> / the raw 2-bit layout keeps
|
||
it O(1).</li>
|
||
<li><strong>Lone-k-mer minimiser.</strong> Routing a <em>synthetic</em> variant to its partition
|
||
needs its minimiser, but <code>RollingStat</code> (<code>obiskbuilder/src/rolling_stat.rs</code>)
|
||
only computes minimisers incrementally along a sequence. Add a standalone
|
||
<code>minimizer(kmer) -> Minimizer</code> that scans the <code>k-m+1</code> m-mer windows
|
||
(<code>PackedSeq::mmer</code>, <code>obikseq/src/packed_seq.rs</code>), canonicalises each, and
|
||
takes the min by <code>seq_hash()</code> — the same selection <code>RollingStat</code> performs,
|
||
evaluated once. Partition index is then
|
||
<code>(minimizer.seq_hash() & (n_partitions - 1)) as usize</code>, exactly as
|
||
<code>QueryBatch::from_records</code> (<code>obikmer/src/cmd/query.rs:142</code>); <code>n_partitions</code>
|
||
is a power of two so the mask is valid.</li>
|
||
</ol>
|
||
<h3 id="step-1-the-tally-accumulator-obikindex">Step 1 — the tally accumulator (<code>obikindex</code>)</h3>
|
||
<p>A <code>SnpTally</code> holding, per genome pair, the 4x4 joint count of central bases:
|
||
<code>n * n * 4 * 4</code> <code>u64</code> (or a packed lower-triangular form since it is
|
||
symmetric). Provide <code>merge(&mut self, other: &SnpTally)</code> for the thread-local
|
||
reduce, and accessors yielding, per pair <code>(i,j)</code>: total off-diagonal (SNP),
|
||
diagonal (shared, i.e. <code>p_hat</code>'s denominator minus SNP), transition count
|
||
<code>P</code>, transversion count <code>Q</code>. The diagonal is always populated — it is not an
|
||
optional LogDet-only extra, since <code>p_hat</code>'s denominator is no longer sourced
|
||
from the external <code>shared_kmers</code> matrix (see "Locus eligibility" and
|
||
"Statistic and correspondence with <code>shared</code>" above): the source k-mer's own
|
||
presence/count vector, already in hand when it is enumerated, supplies the
|
||
diagonal entry directly, at no extra lookup cost.</p>
|
||
<h3 id="step-2-the-sweep-obikindex-new-snprs">Step 2 — the sweep (<code>obikindex</code>, new <code>snp.rs</code>)</h3>
|
||
<p>Mirror <code>distance.rs</code>: open the presence or count store per partition. But
|
||
instead of a per-partition <code>partial_*</code>, run the sequential source sweep:</p>
|
||
<div class="highlight"><pre><span></span><code>for p in 0..n_partitions: # OUTER — sequential
|
||
open source partition p's layers (QueryLayer-style, obikpartitionner)
|
||
enumerate distinct canonical k-mers of p (one per MPHF slot) with their
|
||
presence/count vectors # column-major, as query stage 2
|
||
par_iter over these source k-mers: # INNER — rayon, thread-local tally
|
||
apply eligibility rule to the source's own vector (raw: none; # diagonal
|
||
stringent: exactly one of the 4 forms present in each genome) # gate
|
||
for i in genomes eligible with source base a:
|
||
for j in genomes eligible with source base a:
|
||
thread_tally[i,j][a,a] += 1 # diagonal — no extra lookup
|
||
for each of the 3 central variants:
|
||
q = partition_of(variant)
|
||
if q < p: continue # dedup: forward targets only
|
||
if q == p and variant <= source.raw(): continue # in-partition tie-break
|
||
slot = layers[q].find_slot(variant) # MphfLayer::find, mmap'd
|
||
if hit:
|
||
vb = variant presence/count vector
|
||
apply eligibility rule to vb (as above)
|
||
for i in eligible genomes with source base a:
|
||
for j in eligible genomes with variant base b:
|
||
thread_tally[i,j][a,b] += 1
|
||
merge thread-local tallies into global SnpTally
|
||
</code></pre></div>
|
||
<p>The inner lookup is precisely <code>QueryLayer::find_slot</code> +
|
||
<code>col_value(g, slot)</code> (<code>obikpartitionner/src/query_layer.rs</code>) — reuse or factor
|
||
out that path rather than reimplementing MPHF access. Enumerating "all distinct
|
||
k-mers of a partition with their vectors" is the <code>dump</code>/<code>query</code> stage-2
|
||
column-major scan already implemented in <code>dump_layer.rs</code> /
|
||
<code>query_partition_with</code>; factor a reusable iterator if none fits.</p>
|
||
<p><code>presence_threshold</code> applies exactly as elsewhere: a genome "carries base b"
|
||
iff its count at that slot is <code>>= presence_threshold</code> (trivially <code>>= 1</code> for
|
||
presence indexes).</p>
|
||
<h3 id="open-problem-unresolved-session-end-not-yet-fully-convinced">Open problem (unresolved, session end — not yet fully convinced)</h3>
|
||
<p>The <code>q >= p</code> / tie-break dedup rule in Step 2's pseudocode above is <strong>flawed</strong>
|
||
for the stringent (paralogy-filtered) eligibility rule: it only ever brings
|
||
two family members into view at once (the source and one looked-up variant),
|
||
never all four simultaneously, and which subset gets compared depends on
|
||
partition sweep order. "Exactly one of the 4 forms present in genome A" is a
|
||
whole-family property and cannot be decided correctly from a sequence of
|
||
pairwise, order-dependent glimpses — the pseudocode above needs revision, not
|
||
just the eligibility gate bolted onto it as written.</p>
|
||
<p>Direction discussed, <strong>not yet settled</strong>:</p>
|
||
<ol>
|
||
<li><strong>Every distinct source k-mer looks up all 3 variants unconditionally</strong>
|
||
(drop the <code>q < p</code> skip entirely) so that every observed family member
|
||
independently gathers all 4 vectors (its own + whichever of the 3
|
||
variants exist) at once — a whole-family, order-independent view, computed
|
||
redundantly once per observed member. Same total lookup order of
|
||
magnitude as already budgeted (<code>3 * N_distinct</code>), just organised
|
||
differently (no lookup actually skipped, versus the original rule which
|
||
skipped roughly half).</li>
|
||
<li><strong>Tie-break after gathering, not before</strong>: only the member whose own
|
||
canonical encoding is the smallest <em>among the members actually observed</em>
|
||
(now known, since all were just looked up) writes to <code>SnpTally</code>; the
|
||
others silently discard their redundant computation. Deterministic,
|
||
order-independent — as a side effect this also removes the "outer loop
|
||
must stay sequential" constraint from the cost/parallelism discussion
|
||
above, since no step depends on partition processing order any more.</li>
|
||
<li><strong>Proposed optimisation</strong>: precompute, once at index build time, a
|
||
compact global (not per-genome) annex per MPHF slot — the count of
|
||
<em>other</em> family members observed anywhere in the dataset (0-3). Slots with
|
||
count 0 (majority under low divergence and few genomes, but see the
|
||
scaling caveat below) need no cross-lookup at all: eligibility reduces to
|
||
a local <code>count == 1</code> check at that single slot, and only slots with count<blockquote>
|
||
<p>= 1 enter the 3-lookup sweep machinery above. Revised (see Step 2b
|
||
below): minorant status <em>is</em> stored alongside the count after all, on 3
|
||
bits rather than 2 — it comes for free from the same lookups needed to
|
||
count siblings, and storing it lets the sweep discard non-minorant slots
|
||
without re-fetching anything.</p>
|
||
</blockquote>
|
||
</li>
|
||
</ol>
|
||
<p><strong>Minorant/sibling-count relationship, worked out precisely.</strong> "Minorant" is
|
||
a one-way implication from sibling count, not an equivalence: <code>0 siblings
|
||
=> minorant</code> (trivially — with no other observed member, the k-mer is by
|
||
definition the smallest of the observed set, itself alone), and its
|
||
contrapositive <code>not minorant => >= 1 sibling</code>. The converse does not hold:
|
||
being the minorant says nothing about sibling count — a minorant can have 0,
|
||
1, 2 or 3 siblings, all with larger encodings than itself. Consequence: this
|
||
confirms, as a logical necessity rather than a heuristic, that a 0-sibling
|
||
slot can always write its diagonal contribution with zero ambiguity and no
|
||
lookup (it is unconditionally its own minorant) — but it gives no shortcut
|
||
for the >= 1-sibling case, where minorant status still requires the actual
|
||
comparison of gathered encodings; sibling count alone never determines it.</p>
|
||
<p><strong>When to compute the annex, and cache invalidation.</strong> Sibling count is a
|
||
property of the whole set of columns (genomes/groups) currently in the
|
||
index, not of any single genome — it cannot be computed correctly at
|
||
mono-genome build time (a family may gain siblings, or its minorant may
|
||
change, once more genomes are merged in later). Computing it eagerly at
|
||
every <code>merge</code> would also waste work on intermediate merged states nobody
|
||
ever queries. Instead: compute it lazily, on first <code>phylo</code> call against a
|
||
given index, and persist the result alongside that index for subsequent
|
||
calls — the same lazy-derived-cache pattern <code>PersistentBitMatrix</code> already
|
||
uses for <code>Columnar</code> -> <code>Packed</code>. This requires no explicit invalidation for
|
||
<code>merge</code> or <code>filter</code> (<code>obikindex/src/merge.rs</code>, <code>obikmer/src/cmd/filter.rs</code>):
|
||
both only ever write to a fresh <code>--output</code> directory, never mutate an input
|
||
index in place, so a re-merged/re-filtered index is simply a new state with
|
||
no annex yet. <code>select --in-place</code> (<code>select_layer.rs:139-235</code>) is the
|
||
exception: it aggregates genome columns into groups (Any/All/None/Sum/Min/
|
||
Max) by mutating the existing index's files without changing its location.
|
||
It does not remove k-mer rows, but it can still change eligibility and
|
||
sibling counts derived from those rows (e.g. a <code>Sum</code> over several
|
||
single-copy genomes can read as multi-copy at the group level). Because it
|
||
mutates in place, <strong><code>select --in-place</code> must explicitly invalidate (delete
|
||
or mark stale) any cached sibling-count annex for that index</strong> — the one
|
||
operation in the current pipeline where this doesn't happen for free.</p>
|
||
<p>Not yet convinced this is the right shape, and Step 2's pseudocode above has
|
||
not been rewritten to match — flagged for the next pass rather than resolved
|
||
here.</p>
|
||
<h3 id="step-2b-sibling-count-minorant-annex-consolidated-plan">Step 2b — sibling-count / minorant annex (consolidated plan)</h3>
|
||
<p>Scope: only the precursor annex — not the SNP tally itself, whose Step 2
|
||
sweep remains unresolved above. This piece is simpler than the sweep,
|
||
because it writes to an independent per-slot value, not a shared
|
||
cross-k-mer accumulator, so it needs no dedup/ownership logic at all at this
|
||
stage.</p>
|
||
<p><strong>Revised annex encoding — 4-bit presence mask, not 3-bit (minorant +
|
||
count).</strong> Superseded after settling the "canonical form of a family"
|
||
definition above. The 3-bit design (1 minorant bit + 2-bit sibling count,
|
||
§ below, kept for the historical record) had two problems: it discards
|
||
<em>which</em> variants are present (only how many), so any future consumer
|
||
(the SNP sweep, or a stats pass — see below) that needs to know which bases
|
||
exist still has to regenerate and blindly re-query all 3 candidates; and
|
||
the minorant bit's meaning was tied to whichever member was visited, not to
|
||
a fixed reference. Storing instead a <strong>4-bit mask</strong> — one bit per base
|
||
(A/C/G/T), set iff that member of the family (labelled relative to the
|
||
family's fixed canonical form, i.e. the member with <code>A</code> at the centre — see
|
||
above) is observed anywhere in the index — fixes both:
|
||
- <strong>Sibling count is derived, not stored</strong>: <code>siblings = popcount(mask) - 1</code>.
|
||
- ~~<strong>Minorant is derived, not stored</strong>: regenerate the family's 4 canonical
|
||
forms from the slot's own k-mer (cheap, no lookup — see above), compare
|
||
the raw encodings of whichever bits are set in the mask, take the
|
||
smallest.~~ <strong>Erratum (2026-08-14) — this was wrong, kept struck through
|
||
rather than deleted.</strong> "Cheap, no lookup" only accounts for the bit
|
||
algebra (regenerate 4 forms, compare raw encodings) — true in isolation,
|
||
but it silently assumed "the slot's own k-mer" is a free fact. It isn't:
|
||
getting from an MPHF slot index back to the actual k-mer sequence means
|
||
reconstructing <code>slot_kmer</code> for the whole layer — scan <code>unitigs.bin</code>,
|
||
<code>mphf.find()</code> every k-mer to place it — an O(distinct k-mers) pass
|
||
through the MPHF, not a per-slot O(1) lookup. That reconstruction is free
|
||
<em>only</em> when the caller already needs k-mer identity for something else in
|
||
the same traversal (e.g. the SNP sweep below, which needs it anyway to
|
||
generate <code>central_canonical_neighbors()</code>). A caller that wants <em>only</em> the
|
||
minorant flag pays the full reconstruction for nothing: measured on a real
|
||
run, a bare family-size histogram (four buckets, otherwise near-instant)
|
||
spent 71% of wall-clock in <code>MphfLayer::find</code>, all of it solely to answer
|
||
"is this slot the minorant". <strong>Current design: minorant <em>is</em> stored</strong>
|
||
after all — a 5th mask bit (4 presence bits + 1 minorant bit, still fits
|
||
one byte alongside the presence mask below), written once in the
|
||
construction pass where the k-mer is already in hand for other reasons
|
||
(<code>obikindex::siblings::build_layer_sibling_annex</code>), read back for free by
|
||
every later consumer (<code>FamilyMask::is_minorant</code>). This is genuinely a
|
||
return to the superseded 3-bit design's core idea (store minorant
|
||
alongside the count) — the "which variant" blindness that motivated
|
||
moving away from it is fixed by keeping the full 4-bit presence mask
|
||
too, not by dropping the stored minorant bit again.
|
||
- <strong>A future consumer knows exactly which variants to (re-)query</strong> —
|
||
<code>popcount(mask) - 1</code> lookups instead of always 3, and it knows <em>which</em>
|
||
3 (or fewer) to issue, not just how many hits to expect.
|
||
- The all-zero value (no base present at all) is still logically
|
||
unreachable as a real result — the slot's <em>own</em> base is always present in
|
||
its own family — so it remains available as a free "not yet computed"
|
||
sentinel, exactly as before.</p>
|
||
<ol>
|
||
<li><strong>Primitive.</strong> Reuse <code>central_canonical_neighbors()</code> from Step 0
|
||
unchanged — the 3 canonicalised central-substitution variants of a k-mer
|
||
(plus the identity, i.e. all 4 members of the family — see "Definitions"
|
||
above).</li>
|
||
<li>
|
||
<p><strong>New annex type</strong> (<code>obicompactvec</code>, alongside <code>bitmatrix.rs</code>): a 4-bit-
|
||
per-slot packed array (the presence mask above), one per partition — same
|
||
on-disk shape family as <code>PersistentBitMatrix</code>'s <code>Packed</code> variant, but
|
||
simpler (no per-genome columns, a single derived read-only value per
|
||
slot). As actually implemented (<code>obicompactvec::siblingannex</code>): <strong>not</strong>
|
||
truly bit-packed — 1 byte/slot, 4 presence bits + the minorant bit (see
|
||
the erratum above) in the low 5 bits, 3 unused. Deliberately simpler for
|
||
a first implementation; packing to 5 bits/slot is a pure storage-density
|
||
follow-up, not a behavioural change, still not done as of this note.
|
||
<br />
|
||
<details><summary>Superseded 3-bit design (historical)</summary>
|
||
3 bits, storing minorant status alongside sibling count directly, since
|
||
it came for free from the same lookups (point 3 below) — 5 real states
|
||
(not-minorant; minorant with 0/1/2/3 siblings) fit in 3 bits (8 states,
|
||
3 unused). This let the future SNP sweep discard a non-minorant slot
|
||
instantly, with no lookup at all. The otherwise-unreachable combination
|
||
"not-minorant + 0 siblings" (0 siblings always implies minorant) doubled
|
||
as the "not yet computed" sentinel. Replaced by the 4-bit mask above,
|
||
which subsumes this benefit (minorant still derivable, now for free at
|
||
read time rather than stored) while also fixing the "which variant"
|
||
blindness.
|
||
</details></p>
|
||
</li>
|
||
<li>
|
||
<p><strong>Computation pass</strong> (<code>obikindex</code>, new <code>siblings.rs</code>): <strong>one
|
||
<code>obipipeline</code> run per layer, iterated sequentially over the index's
|
||
layers</strong> — settled after two false starts, worth recording both.</p>
|
||
</li>
|
||
<li><em>False start 1</em>: "fully parallel over every partition/slot at once,
|
||
no ordering at all". Correctness is fine with this (sibling count and
|
||
minorant are order-independent, unlike the old <code>q >= p</code> dedup they
|
||
replace), but it reintroduces, at a larger scale, exactly the
|
||
memory-blowup the original Step 2 sweep's sequential-outer-loop
|
||
constraint existed to prevent: scattering every source partition at
|
||
once multiplies the in-flight outgoing-query volume by the number of
|
||
partitions.</li>
|
||
<li><em>False start 2</em>: push the layer loop itself into the pipeline (source
|
||
= the index's layers, a first <code>Flat</code> stage expands each layer into
|
||
its k-mers). <code>obipipeline</code>'s scheduler already bounds memory on its
|
||
own — it dispatches every item through a <strong>shared</strong> worker pool at
|
||
each stage boundary (<code>scheduler.rs:217-372</code>, <code>dispatch()</code> into a
|
||
common <code>worker_tx</code> queue, any free worker picks up any pending item;
|
||
not "one worker owns a chunk end to end"), with a biased <code>Select</code>
|
||
that prioritises draining items already advanced in the chain over
|
||
admitting new source items (<code>scheduler.rs:271-282</code>: stage results
|
||
outrank the source, "vider le pipeline en priorité" / "dernier
|
||
recours" for new data) — so bounded channel <code>capacity</code> plus this
|
||
drain-first bias already caps in-flight work without any external
|
||
sequential discipline. Correct, but it means k-mers from several
|
||
layers can be completing concurrently, so the sink would need to
|
||
track several open per-layer annex-file writers at once — real,
|
||
avoidable complexity.</li>
|
||
<li><strong>Settled design</strong>: keep the layer loop external and sequential —
|
||
not for memory (the pipeline's own <code>capacity</code>/priority mechanism
|
||
already provides that, for free, regardless), but so each pipeline
|
||
run's sink targets exactly one layer's annex file, no concurrent
|
||
multi-writer bookkeeping. Per layer: source = that layer's distinct
|
||
k-mers; a <code>Flat</code> (1->N) stage generates the 3 central variants of a
|
||
k-mer, each tagged with its origin (local slot); a transform stage
|
||
routes each variant to its target partition (unchanged per-k-mer
|
||
minimiser); a transform stage performs the lookup (existence-only —
|
||
<code>find_slot</code> hit/miss, cheaper than the SNP sweep's full column
|
||
fetch); a final stage/sink folds each answer into its origin's
|
||
running state (below) and, once a layer's k-mers are all resolved,
|
||
flushes the completed array to that layer's annex file. Many small,
|
||
single-purpose stages on purpose, to let the scheduler interleave
|
||
them finely across many in-flight items — this deliberately does
|
||
<strong>not</strong> mirror how <code>obipipeline</code> is used elsewhere today: <code>query.rs</code>'s
|
||
<code>process_chunk</code> lumps parse+route+query+serialise into one closure
|
||
(<code>query.rs:325,743-758</code>), and <code>scatter.rs</code> only pipelines file-
|
||
reading/superkmer construction, routing partitions afterwards in a
|
||
plain sequential loop (<code>KmerPartition::write_batch</code>,
|
||
<code>partition.rs:140</code>) — both under-use the fine-grained scheduling the
|
||
mechanism offers, so they are not precedents to copy, only existing
|
||
(and arguably improvable, out of scope here) usages. Cross-partition
|
||
lookups (querying another layer's MPHF for a variant) remain
|
||
necessary as before — only the <em>output</em> side is kept single-layer.</li>
|
||
<li><strong>Reconciliation</strong>: processed at the granularity of one <em>answer batch
|
||
per destination partition</em>, not one source k-mer at a time — this is a
|
||
proper shuffle, not a per-k-mer wait. Each source partition <code>p</code> holds a
|
||
small array of running states <code>(minorant = true, siblings = 0)</code>, one
|
||
per local slot, initialised at scatter time and <strong>persisting across
|
||
however many destination-partition batches answer it</strong> (up to 3, one
|
||
per variant, not necessarily all from the same <code>q</code>). Every scattered
|
||
query carries an origin tag (source partition + local slot) so its
|
||
answer can be routed back. When target partition <code>q</code> returns its batch
|
||
(all answers for every query that named <code>q</code>, regardless of which source
|
||
k-mer or which source partition they came from), that batch is walked
|
||
once, locally, and each answer updates — via its origin tag — the
|
||
matching entry in <em>its</em> source partition's array: a miss changes
|
||
nothing; a hit does <code>siblings += 1</code>, and if the found sibling's own
|
||
encoding is smaller than the source's, <code>minorant = false</code>. A given
|
||
source k-mer's state is final only once every destination batch
|
||
concerning it has been folded in; its partition's array is flushed to
|
||
the persistent annex once complete. Commutative per entry, so the order
|
||
in which destination batches arrive and get folded in doesn't matter.</li>
|
||
</ol>
|
||
<p><strong>Open optimisation, not adopted yet — real tradeoff, not a free win.</strong>
|
||
Since looking up sibling <code>y</code> from <code>x</code>'s visit already yields everything
|
||
needed to fill <code>y</code>'s own annex entry too, one visit per <em>family</em> could in
|
||
principle replace one visit per <em>observed family member</em> — cutting this
|
||
pass's cost roughly by the average family size instead of paying
|
||
<code>3 * N_distinct</code> regardless. But it means threads processing different
|
||
source k-mers can end up writing the <em>same</em> sibling's slot concurrently —
|
||
the fully independent, ownership-free parallelism of the plan above is
|
||
deliberately traded away for this gain. It stays safe only because the
|
||
computed value for a given slot is deterministic regardless of who
|
||
computes it, so redundant concurrent writes converge to the same
|
||
value — correct as long as each write is atomic, no locking needed — but
|
||
it is a real design complexity increase over "every member redoes its
|
||
own 3 lookups independently," not a strict improvement to adopt by
|
||
default.
|
||
4. <strong>Trigger and caching</strong> (<code>obikindex::KmerIndex</code>/<code>distance.rs</code>): compute
|
||
lazily on first <code>phylo</code> call for an SNP-family metric against a given
|
||
index; check for an existing annex file first (mirrors
|
||
<code>PersistentBitMatrix::open()</code>'s auto-detect-and-fall-back,
|
||
<code>bitmatrix.rs:264-287</code>); if absent, run step 3 and persist; if present,
|
||
mmap and reuse.
|
||
5. <strong>Invalidation.</strong> <code>merge</code> and <code>filter</code> always write to a fresh <code>--output</code>
|
||
directory (<code>obikindex/src/merge.rs</code>, <code>obikmer/src/cmd/filter.rs</code>) so a
|
||
re-merged/re-filtered index simply has no annex yet — nothing to
|
||
invalidate. <code>select --in-place</code> (<code>select_layer.rs:139-235</code>) mutates
|
||
columns of an existing index without changing its location, which can
|
||
change sibling counts without removing rows — it must explicitly delete
|
||
any cached annex for that index as part of its in-place rewrite.
|
||
6. <strong>Testing</strong>: hand-built tiny indexes with known sibling counts (0-3);
|
||
order-independence (recompute twice on a static index, identical
|
||
result, given the fully-parallel no-ownership design); invalidation
|
||
(annex absent/correctly recomputed after <code>select --in-place</code>); once
|
||
Step 2's sweep is fixed, a regression check that sibling_count == 0
|
||
slots are never looked up cross-partition during the sweep.</p>
|
||
<p>Cost: <code>3 * N_distinct</code> existence-only lookups, computed once per index
|
||
state and amortised over every subsequent <code>phylo</code> call that reuses the
|
||
cached annex — cheaper per-lookup than the sweep itself (hit/miss only, no
|
||
column fetch).</p>
|
||
<h3 id="step-3-finalisation-obikindex">Step 3 — finalisation (<code>obikindex</code>)</h3>
|
||
<p>From the global <code>SnpTally</code> alone (diagonal and off-diagonal both populated by
|
||
the sweep, see Step 1/2 — no dependency on the external <code>shared_kmers</code>
|
||
matrix), derive n x n distance matrices, each a pure function of the
|
||
accumulated counts (same shape as <code>jaccard_to_mash</code>):</p>
|
||
<ul>
|
||
<li><code>p_hat[i,j] = SNP / (SNP + shared)</code></li>
|
||
<li>Jukes-Cantor, Kimura-2P (from <code>P</code>, <code>Q</code>), optionally LogDet (needs the
|
||
diagonal + base-composition margins).</li>
|
||
</ul>
|
||
<p>Guard the singularities (<code>p >= 3/4</code> for JC, <code>1-2P-Q <= 0</code> or <code>1-2Q <= 0</code> for
|
||
K2P) by clamping to a max distance, as <code>jaccard_to_mash</code> clamps <code>J <= 0</code>.</p>
|
||
<h3 id="step-4-surfacing-obikindex-obikmer-cli">Step 4 — surfacing (<code>obikindex</code> + <code>obikmer</code> CLI)</h3>
|
||
<p>These metrics do not fit <code>DistanceMetric</code>'s current <code>LayeredStore</code>-partial
|
||
dispatch (they need the cross-partition sweep and produce a different
|
||
intermediate). Two options, to decide:</p>
|
||
<ul>
|
||
<li><strong>(a)</strong> New <code>DistanceMetric</code> variants (<code>Pdistance</code>, <code>JukesCantor</code>,
|
||
<code>Kimura2P</code>, <code>LogDet</code>) whose <code>KmerIndex::distance</code> arm calls the sweep
|
||
(<code>snp.rs</code>) instead of the partial path, still returning <code>DistanceOutput</code>.
|
||
Keeps one CLI surface (<code>--metric jukes-cantor</code>), at the cost of a branch in
|
||
<code>distance()</code> that ignores the <code>LayeredStore</code> it built.</li>
|
||
<li><strong>(b)</strong> A dedicated pathway (<code>KmerIndex::snp_distance</code>) and a distinct CLI
|
||
entry, if mixing a cross-partition sweep into the partition-local <code>phylo</code>
|
||
command is judged architecturally muddy.</li>
|
||
</ul>
|
||
<p>Recommendation: (a) for user ergonomics (all pairwise distances under
|
||
<code>phylo</code>, all feeding NJ/UPGMA/<code>--shared-kmers</code> unchanged), but compute the
|
||
sweep lazily only when an SNP-family metric is requested, so the existing
|
||
metrics keep their partition-local fast path untouched.</p>
|
||
<h3 id="step-5-subsampling-flag">Step 5 — subsampling flag</h3>
|
||
<p>Add <code>--snp-sample <fraction></code> (or a bottom-<code>s</code> hash threshold): restrict the
|
||
source-k-mer enumeration in Step 2 to <code>seq_hash(kmer) < threshold</code>. Divides
|
||
lookups proportionally; <code>p_hat</code> is unbiased. Off by default (exact).</p>
|
||
<h3 id="testing">Testing</h3>
|
||
<ul>
|
||
<li><strong>Primitive unit tests</strong>: <code>central_canonical_neighbors</code> on hand-checked
|
||
k-mers incl. palindrome-boundary cases; lone-k-mer <code>minimizer</code> against
|
||
<code>RollingStat</code>'s incremental result on the same k-mer.</li>
|
||
<li><strong>End-to-end tiny index</strong>: two 1-genome indexes differing by a handful of
|
||
known isolated SNPs (transitions and transversions placed by hand), assert
|
||
exact <code>SNP</code>, <code>P</code>, <code>Q</code> counts and the resulting JC/K2P values.</li>
|
||
<li><strong>Dedup invariant</strong>: assert the tally is identical regardless of genome/
|
||
partition order and that no pair is double-counted (compare against a
|
||
brute-force all-pairs reference on a small index).</li>
|
||
<li><strong>Subsampling</strong>: <code>p_hat</code> within sampling error of the exact run.</li>
|
||
</ul>
|
||
<h3 id="suggested-phasing">Suggested phasing</h3>
|
||
<ol>
|
||
<li>Step 0 primitives + their unit tests (self-contained, no distance wiring).
|
||
This also unblocks the long-declared-but-unimplemented <code>query --mismatch</code>
|
||
(<code>obikmer/src/cmd/query.rs:676</code>, currently a warning), which needs the same
|
||
neighbour + routing machinery.</li>
|
||
<li><code>SnpTally</code> + finalisation math with a brute-force (non-swept) reference
|
||
backend, validated on a tiny index.</li>
|
||
<li>The real per-partition sweep (Step 2) behind the same finalisation; assert
|
||
it matches the brute-force backend.</li>
|
||
<li>CLI surfacing (Step 4a) and NJ/UPGMA integration (already generic over the
|
||
matrix).</li>
|
||
<li>Subsampling (Step 5).</li>
|
||
</ol>
|
||
<h2 id="locus-dropout-under-incomplete-coverage-absence-as-missing-data-not-a-penalized-state-2026-08-13-discussion">Locus dropout under incomplete coverage: absence as missing data, not a penalized state (2026-08-13, discussion)</h2>
|
||
<p>Distinct from the SNP/Sankoff calibration above (which conditions on
|
||
<em>conserved-flank</em> comparability), this addresses the coarser presence/absence
|
||
signal of a family itself, for genome-skim/reduced-representation inputs
|
||
where per-genome coverage is often <code><1x</code>, sometimes <code><0.5x</code>. At that depth,
|
||
non-detection of a locus is dominated by Poisson sampling failure, not by
|
||
true loss — even for a single-copy locus, <code>P(non-detection | present)</code> is
|
||
non-negligible. Treating "absent" as a real 0-state therefore risks
|
||
clustering genomes by shared undersampling rather than shared ancestry.</p>
|
||
<p><strong>Parsimony.</strong> Loss and gain are the two directions of the same transition,
|
||
so they cannot be freed asymmetrically without a biological justification for
|
||
trusting one direction and not the other — none exists here. Freeing both
|
||
directions to zero cost makes the character contribute nothing to any tree's
|
||
score regardless of topology, which is operationally equivalent to just
|
||
coding the state as ambiguous (<code>?</code>) rather than 0/1. <code>?</code> is the more honest
|
||
encoding when downstream ancestral-state reconstruction is wanted: cost-zero
|
||
Fitch/Sankoff still assigns internal states arbitrarily among equal-cost
|
||
options, whereas <code>?</code> makes the "we don't know" explicit.</p>
|
||
<p><strong>ML.</strong> The correct analogue is not a tuned asymmetric loss/gain rate but
|
||
literal missing-data coding at the tip: an ambiguous state (all-states tip
|
||
likelihood vector, e.g. <code>(1,1)</code> for a 2-state character) so Felsenstein
|
||
pruning marginalises over it and the character contributes likelihood 1
|
||
(neutral) for that taxon — the same mechanism already used for the <code>N</code>
|
||
IUPAC-ambiguous case in the 16-state Sankoff encoding (see "IQ-TREE
|
||
custom-model format" above). No rate parameter to justify or estimate.</p>
|
||
<p><strong>Open item, not yet decided.</strong> No per-locus/per-genome confidence
|
||
distinction is implemented: the current plan is to code <em>every</em> non-detected
|
||
locus as missing/ambiguous, uniformly, without a depth-based threshold for
|
||
"confidently absent." Deferred because at typical skim coverage the
|
||
confidently-absent case is rare enough that a uniform treatment is
|
||
defensible as a first pass. Known consequence to watch for: a genome with
|
||
very low overall coverage ends up with a mostly-<code>?</code> row and constrains the
|
||
tree weakly — will likely need a minimum-detected-loci-per-genome filter
|
||
upstream (exclude or flag under-covered genomes) before this becomes a
|
||
correctness problem rather than a power problem. Not designed yet.</p>
|
||
<h3 id="-free-loss-implemented-2026-08-13"><code>--free-loss</code> implemented (2026-08-13)</h3>
|
||
<p><code>obikmer phylo --sankoff --free-loss</code> (with <code>--tnt</code>/<code>--phyg</code>/<code>--iqtree</code>)
|
||
recodes <code>∅</code> (a family's non-detection in a genome) to <code>?</code> instead of the
|
||
default ordinary, costed 16th state — the resolution of the discussion
|
||
above. <code>?</code>, not <code>-</code>: <code>-</code> still carries gap/indel semantics in TNT/PhyG/
|
||
IQ-TREE (see "<code>∅</code> stays an ordinary 16th state" above, and the reasoning
|
||
that motivated it in the first place — the RAxML-era failure where <code>-</code>
|
||
was silently swallowed as missing data was undesired back when <code>∅</code> was
|
||
meant to be fully informative; <code>--free-loss</code> is the first case where that
|
||
same tool behaviour is exactly what's wanted). <code>?</code> is the standard
|
||
missing/ambiguous symbol these tools already recognise without needing to
|
||
be declared in the custom alphabet, so parsimony treats it as
|
||
"compatible with any state" (Sankoff/Fitch's own missing-data extension —
|
||
no per-transition cost edit needed) and Felsenstein pruning marginalises
|
||
it for IQ-TREE's ML.</p>
|
||
<p>Implementation is a recoding at each export's existing <code>-</code>-normalisation
|
||
boundary, not a change to matrix calibration: <code>write_sankoff_alignment_fasta</code>
|
||
(<code>sankoff.rs</code>, feeds <code>--phyg</code> too via its <code>prefasta:</code>), <code>write_sankoff_tnt</code>
|
||
(<code>tnt.rs</code>, <code>?</code> written directly, bypassing <code>TNT_STATE_SYMBOL</code> — there is no
|
||
state index for it), and <code>compact_alphabet</code>/<code>write_iqtree_alignment</code>
|
||
(<code>iqtree.rs</code>, <code>∅</code> excluded from the compact alphabet entirely rather than
|
||
assigned a symbol). The calibrated cost matrix itself is unchanged — its
|
||
<code>∅</code>-row/column simply becomes unused once no sequence contains that state.</p>
|
||
<p><strong>Unverified, flagged rather than assumed</strong> (this project's own standard for
|
||
tool-specific claims — see the TNT/PhyG/IQ-TREE syntax facts above, all
|
||
checked against the real binaries, not the manuals alone):
|
||
- PhyG auto-adds its own <code>(n+1)</code>-th gap/indel row to the <code>tcm:</code> matrix
|
||
(<code>phyg.rs</code>); with <code>--free-loss</code> the reused <code>_sankoff.fasta</code> now contains
|
||
real <code>?</code> characters PhyG did not have to interpret before. Not confirmed
|
||
against the local PhyG binary that <code>?</code> is read as "any state" rather than
|
||
routed through that gap row, or otherwise mishandled.
|
||
- IQ-TREE's exact <code>?</code>/missing-data semantics for <code>--seqtype MORPH</code> with a
|
||
custom-file model plus <code>+ASC</code> — assumed to marginalise via Felsenstein
|
||
pruning by the general convention, not yet re-verified against the local
|
||
<code>iqtree3</code> build the way the rest of the custom-model mechanism was.</p>
|
||
<p>Both should be checked against a real run (small index) before <code>--free-loss</code>
|
||
results are trusted for TNT/PhyG/IQ-TREE, the same way every other
|
||
tool-specific behaviour in this file was established.</p>
|
||
<h3 id="tnt-confirmed-2026-08-13-iq-tree-asc-constant-column-fix-2026-08-13">TNT confirmed (2026-08-13); IQ-TREE <code>+ASC</code> constant-column fix (2026-08-13)</h3>
|
||
<p><strong>TNT.</strong> Real run, 20-genome benchmark index, <code>--free-loss --tnt</code>: <code>?</code>
|
||
present for all 20 taxa in both <code>sankoff.fasta</code> and <code>sankoff.tnt</code>, loaded
|
||
with no syntax error and no triangle-inequality warning, <code>mult</code> search
|
||
completed, <code>sankoff.tre</code> holds 20 equally-parsimonious trees. Topology is
|
||
exactly what the theory predicts: strain-level clades unaffected (the 4 <em>E.
|
||
coli</em>, the <em>Klebsiella</em> trio, the <em>Salmonella</em> pairs — driven by real
|
||
substitution signal, not gene-content) but resolution among the deeply
|
||
divergent lineages (<em>Candidozyma</em>, <em>Saccharolobus</em>, <em>Wolbachia</em>,
|
||
<em>Acidobacterium</em>, <em>Opitutus</em>, <em>Bacillus</em>/<em>Shouchella</em>) is weak and unstable
|
||
across the 20 tied trees — the expected trade-off of discarding the
|
||
gene-content/dropout signal that used to help place exactly those lineages.
|
||
TNT's own handling of <code>?</code> — read as "compatible with any state" — no longer
|
||
a flagged assumption.</p>
|
||
<p><strong>IQ-TREE: <code>+ASC</code> broke on real data — non-informative columns survive
|
||
<code>--free-loss</code>, distinct bug from the <code>--exclude-genome</code> one above.</strong>
|
||
<code>snp_pseudo_alignment</code>'s "variable family" test (<code>family_size() >= 2</code>) is a
|
||
whole-annex property, blind to any single column's actual calls. Once
|
||
non-detection is hidden as <code>?</code> (ignored for constancy) rather than counted
|
||
as its own real state, a family that is genuinely variable <em>somewhere in the
|
||
dataset</em> can still have only one distinct call among the genomes that
|
||
actually carry it in one particular column, with the rest <code>?</code> — constant
|
||
once missing data is excluded, which is exactly what <code>+ASC</code> forbids. Same
|
||
underlying failure mode as <code>--exclude-genome</code>'s (documented above, "Two
|
||
consistency bugs found and fixed post-implementation"), same remedy: rescan
|
||
columns and drop the ones that are constant once <code>-</code> is ignored — but a
|
||
different trigger (hiding cells vs. dropping whole rows), so a separate fix
|
||
was needed rather than reusing <code>drop_excluded</code> directly. Implemented as
|
||
<code>iqtree::drop_ascertainment_noninformative</code>, run only under <code>--free-loss</code>,
|
||
only on IQ-TREE's own copy of the alignment (<code>--tnt</code>/<code>--phyg</code> have no
|
||
no-invariant-site requirement and must not lose those columns). Not yet
|
||
re-run against the local <code>iqtree3</code> binary to confirm <code>+ASC</code> now accepts the
|
||
filtered alignment — next verification step, alongside the still-open PhyG
|
||
<code>?</code>-vs-gap-row question above.</p>
|
||
<h3 id="-family-overlap-min-shared-family-and-a-second-free_loss-bug-cardinality-transition-cost-not-zeroed-2026-08-13"><code>--family-overlap</code>/<code>--min-shared-family</code>, and a second <code>free_loss</code> bug: cardinality-transition cost not zeroed (2026-08-13)</h3>
|
||
<p><strong>Diagnostic added: <code>--family-overlap</code>.</strong> Real per-genome-pair count of
|
||
shared variable families (columns where neither genome is <code>∅</code>), diagonal
|
||
kept (degenerates to "genome's own total family count", not a separate
|
||
computation — see <code>obikmer/src/cmd/phylo/family_overlap.rs</code>). On the
|
||
20-genome benchmark this predicted the <code>--tnt</code>/<code>--iqtree</code> rogue-placement
|
||
pattern exactly: Candidozyma shares <strong>0</strong> families with all 19 others (hence
|
||
its ~0-length branch, placed anywhere); the mean off-diagonal shared-family
|
||
count per genome splits cleanly below/above ~1000 (<code>Candidozyma</code> 0,
|
||
<code>Saccharolobus</code> 0.9, <code>Opitutus</code> 20.5, <code>Wolbachia</code> 37.1, <code>Acidobacterium</code>
|
||
41.2, <code>Shouchella</code> 49.4, <code>Bacillus</code> 54.5, <code>Proteus</code> 629, <code>Yersinia</code> 907.5 —
|
||
all below 1000 — vs. <code>Klebsiella</code> ~5000-6000, <code>Salmonella</code> ~8000-17000,
|
||
<code>E. coli</code> ~19000-24000, all comfortably above). <code>--min-shared-family N</code>
|
||
automates exclusion below that threshold, reusing the <code>--exclude-genome</code>
|
||
machinery. Restricting to the 3 genera above threshold (11 genomes: 4
|
||
<em>E. coli</em>, 3 <em>Klebsiella</em>, 4 <em>Salmonella</em>) gives 3 cleanly separated
|
||
monophyletic clades in both TNT and IQ-TREE (without <code>+ASC</code>, see below) —
|
||
confirms the family-overlap statistic, not just the tree topology, as the
|
||
right lens for this failure mode. <code>N=1000</code> is benchmark-specific (huge
|
||
cross-domain divergence), not a universal constant — recompute
|
||
<code>--family-overlap</code> and read its own gap on any other dataset.</p>
|
||
<p><strong>Second <code>free_loss</code> bug, found by re-deriving the design rather than by a
|
||
failing run: the cardinality-transition cost wasn't zeroed.</strong>
|
||
<code>--free-loss</code> neutralises <code>∅</code> (whole-family non-detection) by recoding it to
|
||
<code>?</code>, so transitions to/from the empty state are never scored. But the
|
||
calibrated 16-state cost matrix (<code>obikindex::pairwise_cost_matrix</code>) prices
|
||
<em>every</em> cardinality change — not just <code>X → ∅</code> — via a shared
|
||
<code>P_cardinality(|A|→|B|)</code> factor (the same constant this project has called
|
||
<code>c_ctx</code> throughout: "used everywhere a member is gained or lost", see "A
|
||
concrete Sankoff cost matrix" above). So a transition between two
|
||
<em>non-empty</em> states that differ in how many family members are detected
|
||
(e.g. <code>{A} -> {A,C}</code>, gaining a second sibling) was still carrying the full
|
||
calibrated gain/loss penalty under <code>--free-loss</code>, even though the same
|
||
low/incomplete-coverage argument for not trusting whole-family absence
|
||
applies identically one level down: whether a genome shows 1 vs 2 detected
|
||
members of a family it does carry is exactly as vulnerable to sampling
|
||
failure as whether the family was detected at all. Fixed by adding a
|
||
<code>free_loss: bool</code> parameter to <code>pairwise_cost_matrix</code> that drops the
|
||
<code>P_cardinality</code> factor from the log-probability sum entirely when set,
|
||
leaving cost driven only by composition matching (shared-base retention,
|
||
paired substitutions) — never by cardinality difference alone. Unit-tested
|
||
(<code>free_loss_ignores_cardinality_transition_probs</code>): two wildly different
|
||
cardinality models produce an identical cost matrix once <code>free_loss</code> is
|
||
true. Not yet re-run against a real index to see how much this changes the
|
||
resulting tree beyond the unit-level guarantee.</p>
|
||
<p><strong><code>+ASC</code> misbehaves on closely-related taxa, separately from either bug
|
||
above.</strong> Restricting to the 3-genus, 11-genome subset and running <code>--iqtree</code>
|
||
with <code>+ASC</code> converges suspiciously fast compared to the same run without
|
||
<code>+ASC</code> (which takes the expected ~100 NNI iterations, ~1m38s, converges to a
|
||
stable score, total tree length 7.147, and gives the clean 3-clade
|
||
topology). Not yet root-caused with a real <code>+ASC</code> log in hand (the one run
|
||
overwrote its own log under the shared <code>--prefix</code>, see the operational note
|
||
below) — working hypothesis, not confirmed: <code>+ASC</code>'s correction divides by
|
||
<code>P(site variable | current model, current branch lengths)</code>, recomputed
|
||
every iteration; on closely-related strains (short branches), a plausible
|
||
model predicts most sites <em>should</em> be constant, so <code>P(variable)</code> can be
|
||
small enough that the correction becomes numerically unstable. Needs a
|
||
dedicated run (own <code>--prefix</code>, log kept) before this is more than a
|
||
hypothesis.</p>
|
||
<p><strong>Operational note, not a code bug: IQ-TREE's <code>--undo</code>/checkpoint reuse
|
||
across different alignments crashes.</strong> Reusing <code>iqtree.ckp.gz</code> (and the tree
|
||
it stores) from a run on one taxon set against a regenerated <code>iqtree.fasta</code>
|
||
with a different taxon count crashes with <code>Assertion 'i >= 0 && i <
|
||
(int)seq_names.size()' failed</code> in <code>Alignment::getSeqName</code> /
|
||
<code>PhyloTree::assignLeafNames</code> — <code>--undo</code> tries to graft a stored tree with
|
||
the old leaf count onto the new, smaller alignment. Not an <code>obikmer</code> bug;
|
||
avoid by using a fresh <code>--prefix</code> (or <code>-redo</code>) whenever the underlying
|
||
alignment's taxon set changes, never <code>--undo</code> across them.</p>
|
||
<h2 id="references">References</h2>
|
||
<p>The Mash mutation-rate model this discussion contrasts with:
|
||
(Fan <em>et al.</em> 2015; Marbl Lab 2026)<sup id="fnref:Mash-distances-doc"><a class="footnote-ref" href="#fn:Mash-distances-doc">1</a></sup> <sup id="fnref:Fan2015-mash-formula"><a class="footnote-ref" href="#fn:Fan2015-mash-formula">2</a></sup>.</p>
|
||
<div class="footnote">
|
||
<hr />
|
||
<ol>
|
||
<li id="fn:Mash-distances-doc">
|
||
<p>Marbl Lab. (2026). <a href="https://mash.readthedocs.io/en/latest/distances.html">Mash distance</a>. <a class="footnote-backref" href="#fnref:Mash-distances-doc" title="Jump back to footnote 1 in the text">↩</a></p>
|
||
</li>
|
||
<li id="fn:Fan2015-mash-formula">
|
||
<p>Fan, H., Ives, A.R., Surget-Groba, Y. & Cannon, C.H. (2015). <a href="https://doi.org/10.1186/s12864-015-1647-5">An assembly and alignment-free method of phylogeny reconstruction from next-generation sequencing data</a>. <em>BMC Genomics</em>, 16. <a class="footnote-backref" href="#fnref:Fan2015-mash-formula" title="Jump back to footnote 2 in the text">↩</a></p>
|
||
</li>
|
||
</ol>
|
||
</div>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
</article>
|
||
</div>
|
||
|
||
|
||
<script>var target=document.getElementById(location.hash.slice(1));target&&target.name&&(target.checked=target.name.startsWith("__tabbed_"))</script>
|
||
</div>
|
||
|
||
</main>
|
||
|
||
<footer class="md-footer">
|
||
|
||
<div class="md-footer-meta md-typeset">
|
||
<div class="md-footer-meta__inner md-grid">
|
||
<div class="md-copyright">
|
||
|
||
|
||
Made with
|
||
<a href="https://squidfunk.github.io/mkdocs-material/" target="_blank" rel="noopener">
|
||
Material for MkDocs
|
||
</a>
|
||
|
||
</div>
|
||
|
||
</div>
|
||
</div>
|
||
</footer>
|
||
|
||
</div>
|
||
<div class="md-dialog" data-md-component="dialog">
|
||
<div class="md-dialog__inner md-typeset"></div>
|
||
</div>
|
||
|
||
|
||
|
||
|
||
|
||
<script id="__config" type="application/json">{"annotate": null, "base": "../..", "features": [], "search": "../../assets/javascripts/workers/search.2c215733.min.js", "tags": null, "translations": {"clipboard.copied": "Copied to clipboard", "clipboard.copy": "Copy to clipboard", "search.result.more.one": "1 more on this page", "search.result.more.other": "# more on this page", "search.result.none": "No matching documents", "search.result.one": "1 matching document", "search.result.other": "# matching documents", "search.result.placeholder": "Type to start searching", "search.result.term.missing": "Missing", "select.version": "Select version"}, "version": null}</script>
|
||
|
||
|
||
<script src="../../assets/javascripts/bundle.79ae519e.min.js"></script>
|
||
|
||
<script src="https://unpkg.com/mathjax@3/es5/tex-mml-chtml.js"></script>
|
||
|
||
|
||
</body>
|
||
</html> |