Skip to content

Phylogeny

Methods for estimating and operating on phylogenies.

Table of Contents

Shear (subset to tips)

Subset ("shear" / prune) a tree down to a specified set of tips. This is the operation you want when a large reference phylogeny is persisted (e.g. in Parquet via COPY ... (FORMAT NEWICK)) and you need to restrict it to the tips relevant to one analysis.

Function signature:

shear_tree(tree_table, tips_table, collapse := true, ignore_missing := false)

Parameters: - tree_table (VARCHAR): Name of a table or view containing tree data in read_newick schema (node_index, parent_index, name, branch_length, edge_id). - tips_table (VARCHAR): Name of a table or view with a name (VARCHAR) column listing the tip names to keep. Passing the tips as a table (rather than a list literal) keeps the call practical when the kept set is large. Duplicate and NULL names are ignored. - collapse (BOOLEAN, default true): How to treat internal nodes left with a single child after pruning. - true — standard phylogenetic shear: remove single-descendant ancestors and sum their branch lengths onto the surviving descendant edge, so tip-to-tip distances are preserved. The lowest common ancestor (LCA) of the kept tips becomes the new root; nodes above it are dropped. - false — preserve every internal node that lies on a path to a kept tip, with original parent links and branch lengths unchanged (unifurcations retained). - ignore_missing (BOOLEAN, default false): When false, a requested tip name not present as a tip in the tree is an error (the message lists the missing names). When true, absent names are skipped. Either way, if no requested tip matches, the call errors (a tree cannot be sheared to nothing).

Output schema: Same as read_newick (without filepath), reindexed 0-based: - node_index (BIGINT), name (VARCHAR), branch_length (DOUBLE, nullable), edge_id (BIGINT, nullable), parent_index (BIGINT, nullable, NULL for root), is_tip (BOOLEAN).

The schema is UNION ALL-compatible with read_newick and can be written straight back out with COPY ... (FORMAT NEWICK).

Behavior / semantics: - Only tip names are matched — an internal-node label in tips_table is treated as missing. - Branch-length summation (collapse mode) treats an unspecified length (NaN → NULL) as 0, but a collapsed chain that is entirely unspecified stays NULL, so topology-only trees (cladograms) round-trip as cladograms. - A collapsed edge keeps the surviving (lower) node's edge_id; the edge_ids of merged intermediate nodes are dropped. - The new root has no incoming edge, so its branch_length is NULL. - Reads both tree_table and tips_table via a fresh connection (works with tables and views); single-threaded build.

Examples: 

-- Persisted reference tree in Parquet, and the tips this analysis cares about.
CREATE TABLE ref_tree AS SELECT * FROM 'reference_tree.parquet';
CREATE TABLE keep AS SELECT feature_id AS name FROM my_feature_table;

-- Standard shear: collapse single-child ancestors, preserve distances.
SELECT * FROM shear_tree('ref_tree', 'keep');

-- Keep every internal node on the retained paths instead.
SELECT * FROM shear_tree('ref_tree', 'keep', collapse := false);

-- Shear many trees against one shared tip list, tolerating tips absent from a
-- given tree.
SELECT * FROM shear_tree('ref_tree', 'keep', ignore_missing := true);

-- Write the sheared tree back out to Newick.
COPY (
    SELECT node_index, name, branch_length, edge_id, parent_index
    FROM shear_tree('ref_tree', 'keep')
) TO 'sheared.nwk' (FORMAT NEWICK);

Error conditions: - tree_table or tips_table does not exist. - tree_table missing required node_index / parent_index columns; tips_table missing required name column. - A requested tip is not a tip in the tree and ignore_missing := false. - No requested tip matches any tip in the tree.

Performance note (large trees): building and shearing a large tree is allocation-bound — profiling a multi-million-node shear shows ~40% of the time in the system allocator (many small allocations for node names, child lists, and index maps). glibc's allocator is comparatively slow for this pattern. Because miint runs inside a host process, the safe way to switch allocators is a whole-process launch option rather than anything baked into the extension: preload jemalloc so every allocation in the process uses it uniformly.

# Linux: ~30% faster on large tree build/shear in our benchmarks.
LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libjemalloc.so.2 duckdb   # or python, etc.

This applies to any allocation-heavy miint workload (large read_newick builds, tree_resolve_placement, QC), not just shear_tree. It does not reduce peak memory materially — it is a speed optimization.

Resolve multifurcations

Resolve every multifurcation (a node with more than two children — a "polytomy") in a tree into a series of bifurcations, producing a strictly bifurcating tree. This is the operation you want before running a method that requires a binary tree — notably phylo_independent_contrasts — on a tree that contains polytomies (e.g. an unrooted tree whose root is a trifurcation, as phylogeny_fasttree emits, or a tree flattened from a taxonomy).

Function signature:

tree_resolve_multifurcations(tree_table)

Parameters: - tree_table (VARCHAR): Name of a table or view containing tree data in read_newick schema (requires node_index and parent_index; name, branch_length, edge_id optional).

Output schema: Same as read_newick (without filepath), reindexed 0-based: node_index, name, branch_length, edge_id, parent_index, is_tip. The schema is UNION ALL-compatible with read_newick, can be written back out with COPY ... (FORMAT NEWICK), and feeds directly into phylo_independent_contrasts.

Behavior / semantics: - A node with m ≥ 3 children c0, c1, …, c(m-1) (in node_index order) is resolved into a deterministic left-comb: the node keeps c0 and a new internal node N1; N1 keeps c1 and N2; … ; N(m-2) keeps c(m-2) and c(m-1). This inserts m-2 new internal nodes. - New connector nodes are unnamed, have branch length 0, and no edge_id. Because the connectors are zero-length, the original edge lengths and all root-to-tip distances are unchanged. - Nodes with two or fewer children are left untouched. Single-child unifurcations are not collapsed — that is a different operation; use shear_tree with collapse := true for that. - A multifurcating root is resolved as well. - Child order follows node_index (the canonical order read_newick assigns, encoding the original Newick left-to-right order); the reader sorts input rows by node_index, so the result is reproducible regardless of how the tree relation is physically stored.

Important — the resolution is statistically arbitrary. A polytomy carries no information about the order in which its lineages diverged, so any bifurcating resolution is as valid as any other. When the resolved tree feeds a downstream analysis such as independent contrasts, each artificially-created dichotomy contributes a contrast that is not real independent information; the effective degrees of freedom should be reduced accordingly. Resolving is a pragmatic step to satisfy a binary-tree requirement, not a reconstruction of unknown branching order.

Examples: 

-- A tree with a trifurcating (unrooted) root, e.g. from phylogeny_fasttree.
CREATE TABLE ref_tree AS SELECT * FROM read_newick('unrooted.nwk');

-- Resolve polytomies into bifurcations (zero-length connectors).
SELECT * FROM tree_resolve_multifurcations('ref_tree');

-- Prep a tree for PIC: resolve, persist, then run independent contrasts.
CREATE TABLE binary_tree AS SELECT * FROM tree_resolve_multifurcations('ref_tree');
SELECT * FROM phylo_independent_contrasts('binary_tree', 'traits');

-- Write the resolved tree back out.
COPY (
    SELECT node_index, name, branch_length, edge_id, parent_index
    FROM tree_resolve_multifurcations('ref_tree')
) TO 'binary.nwk' (FORMAT NEWICK);

Error conditions: - tree_table does not exist. - tree_table missing required node_index / parent_index columns.

Phylogenetic placement (krepp)

Place query sequences onto a backbone phylogeny using krepp, which matches k-mers against a prebuilt LSH index and maximises a pseudo-likelihood over candidate edges. Column names follow read_jplace where the two overlap, so downstream SQL written for one mostly reads the other — but the schemas are not union-compatible: read_jplace ends in filepath and reports only the best placement per fragment, while this ends in distance and, with multi on, reports every candidate edge.

Requires an index built beforehand, either by krepp_index_create from relations already in the database, or by the krepp command-line tool (krepp index -i <map.tsv> -o <index_dir> -t <tree.nwk>).

Function signature:

place_krepp(query_table, index_path [, newick_path, hdist_th, tau, chisq, multi, filter])

Parameters: - query_table (VARCHAR, required): Name of a table or view holding the query sequences. Requires read_id and sequence1 columns, the same contract as the aligners; read_id may be VARCHAR, BIGINT or UUID. - index_path (VARCHAR, required): Path to a krepp index directory. - newick_path (VARCHAR): Backbone tree in Newick format. Required only when the index was built without a guide tree (krepp index without -t); otherwise the tree stored in the index is used. When the Newick carries jplace-style {N} edge decorations those numbers are reported verbatim, which is what lets placements accumulate against a phylogeny whose edge IDs are assigned elsewhere. krepp requires the decoration on every node or none.

krepp's Newick reader is strict, and reports most violations by calling error_exit — std::exit, which would terminate the DuckDB process rather than raise a SQL error. One it does not report at all: it quietly folds stray whitespace into the following taxon name, and the renamed tip then fails to match the index and is dropped without a word. place_krepp checks the tree itself first and raises in both cases. Trees it rejects, all of which many other tools accept — indented and pretty-printed Newick in particular, which is what several tools emit by default:

Rejected Why it matters
Anything after the final ; beyond one newline A trailing blank line, or CRLF line endings
An unquoted [ or ] A [&R] rooted marker, or an inline [comment]
A ; before the end of the file Concatenated trees, e.g. bootstrap replicates
Unquoted whitespace An indented or line-wrapped tree. Also the quiet one: krepp folds a space into the following name, so (A:1, B:1) yields a tip called " B" that silently fails to match the index
A unifurcation (a node with exactly one child) krepp requires every internal node to have two or more children
Two nodes sharing a name krepp reads an internal node's bootstrap support value as its name, so a backbone carrying support collides wherever two clades share a value. Strip support values first

Two more, both about jplace {N} edge decorations, which krepp also enforces by error_exit: a partially decorated tree (krepp wants all nodes or none), and the same edge number on two nodes. Both are pre-checked as well.

A third is miint's own. An edge number outside [0, 4294967295] is rejected here even though the Newick is well formed, because krepp stores it as a 32-bit unsigned value: {-5} would come back as 4294967291, which would break the "reported verbatim" promise above, and {-1} would collide with {4294967295} inside krepp while looking distinct to the uniqueness check.

The tip labels have to be the index's reference names. This is the one failure krepp does not report at all — not even by exiting. A backbone leaf it cannot find in the index is skipped silently, and every skipped leaf makes its edges ineligible, so a tree whose labels don't match returns zero rows with no error: indistinguishable from "none of these reads placed". place_krepp compares the overlap itself and:

Overlap Behavior
The tree has no labeled tips at all Raises. Nothing can be matched by name
No tip in common with the index Raises. The query could not have returned anything
Some tips in common Warns, naming the counts, and places on the matched part of the tree. A backbone covering more than one index is legitimate
All tips in common Silent

Write the tree on one line, with no comments, no support values, and a single trailing newline.

One known limitation. miint reads a {N} decoration only where the jplace specification puts it — after the branch length, (A:1{0},B:1{1}):0{2};. krepp additionally accepts it straight after the closing parenthesis, ((A:1,B:1){1}:1);. A backbone in that second form is rejected here even though krepp would read it. The fix is to widen miint's Newick parser, which is shared with read_newick and every other phylogeny function, so it is not this function's to make. - hdist_th (UINTEGER, default 4): Maximum Hamming distance for a k-mer to count as a match. Must not exceed the index's k, and is rejected at bind time above 31 (the largest k krepp will build) and at scan start above the index's actual k. krepp itself imposes no upper bound and reads past the end of an internal table when one is exceeded — at hdist_th = 4294967295 that was measured to crash the process. - tau (UINTEGER, default 2): Hamming distance threshold for the placement filter. Must not exceed hdist_th, and is rejected at bind time if it does. krepp's CLI refuses the same combination; the library path this function uses would instead read past the end of a hdist_th + 1 element histogram. - chisq (DOUBLE, default 2.706): Chi-square cutoff for distinguishing candidate edges. Must be positive, and is rejected at bind time if it is not (krepp's CLI applies the same constraint); a non-positive cutoff discards every candidate and returns nothing. The default is the 90th percentile of the chi-square distribution with 1 degree of freedom, i.e. a significance level of 0.10. (krepp's own help text calls this "alpha = 90%"; alpha conventionally names the significance level, so read it as the percentile.) - multi (BOOLEAN, default true): Report every candidate edge rather than only the best one. multi=false is markedly less reproducible than the default — see Reproducibility below. - filter (BOOLEAN, default true): Drop placements without enough k-mer support.

Output schema: - fragment (same type as the query relation's read_id: VARCHAR, BIGINT or UUID): The read_id of the placed query, so it joins back to query_table without a cast - edge_num (BIGINT): Edge of the backbone tree the placement sits on - likelihood (DOUBLE): Log likelihood of the placement - like_weight_ratio (DOUBLE): Likelihood weight ratio. Sums to 1 across one fragment's rows with multi on; with multi=false the single row keeps its share of a total computed over all the candidates that were discarded, so it is below 1 whenever there was more than one. It is a distribution over the candidate edges of one read, so aggregate it per fragment (or per edge_num across reads) — a bare sum(like_weight_ratio) over the whole table is not a meaningful quantity. - distal_length (DOUBLE): Distance from the distal (child) end of the edge to the attachment point, always within the edge. With t the edge's branch length (0 if it has none), d_c the read's Jukes-Cantor distance to the edge's child node, and d_p its distance to the parent node: t·d_c/d_p, clamped strictly inside the edge (or 0 when t is 0), when d_p is finite and positive; otherwise min(t/2, d_c) (krepp's make_placement, ext/krepp/src/query.cpp:218-242). krepp has a finite d_p when it estimated a distance for the parent node, which does not require the parent's own edge to be among the output rows; with a single matched node it has none. So distal_length depends only on t, d_c and d_p: reads on one edge get the same value when those agree, and also whenever they land on the clamp next to the parent end (d_c ≥ d_p), on the fallback with d_c ≥ t/2, or on an edge whose length is 0 or missing. - pendant_length (DOUBLE): Length of the branch joining the read to the attachment point, never negative: max(0, min(d_c − distal, d_p − t + distal)) in the first case above, max(0, d_c − distal) in the second. - distance (DOUBLE): krepp's estimated distance to the placement. Not part of jplace; read_jplace has no equivalent column.

Behavior: - One row per candidate edge, so rows are not one-to-one with reads. With multi off there is at most one row per read. - Reads that krepp cannot support are absent from the output entirely, so fragment is a subset of the input read_ids. - Sequences shorter than the index's k are skipped, with one warning per place_krepp call — not one per skipped sequence (see miint_warnings). A sequence containing anything outside the IUPAC nucleotide alphabet raises rather than being silently dropped. - RNA input is accepted: U/u are rewritten to T/t before krepp sees them, because krepp treats uracil as an ambiguous base and an untouched RNA read would place nothing. - Reads are placed on all of DuckDB's threads, and rows come out in no fixed order. Add an ORDER BY when the order matters. - The index is loaded once per call and shared across that call's threads; indexes routinely run to tens of gigabytes. Two place_krepp calls in one statement load two copies concurrently, so a query that compares several parameter settings side by side multiplies the memory accordingly.

Reproducibility: krepp's output is not bit-reproducible, and neither is this function. krepp accumulates and normalizes placements through hash maps keyed by tree-node pointer, so iteration order moves with the heap layout. What that costs you, measured on the toy index:

The unit that matters is the index load, not the process: each place_krepp call loads its own copy in InitGlobal, so two calls in one session are already in different regimes. Everything below was re-measured on krepp v0.10.0 (31205033) against the toy index that run_tests.sh builds, and every figure is an observation on that index rather than a bound.

krepp v0.9.1 added a de-duplication pass to its k-mer match accumulation, described upstream as a consistent tie break. Measured on v0.9.1, it did not make the output reproducible; the figures were essentially unchanged from v0.9.0.

Quantity Reproducible?
Which reads place at all (fragment set), and the total row count Stable in every run measured — 5 separate processes, identical (404 rows over 92 fragments)
likelihood, distance, distal_length, pendant_length Bit-identical across repeated calls, and across 5 separate processes
Which candidate edges appear (the (fragment, edge_num) set) Identical across repeated calls
like_weight_ratio No. Three processes were each asked for the same placements twice; the two answers differed in 59, 55 and 73 of 404 rows respectively, by up to 2.5e-4 absolute. Per-fragment sums still come to 1 (measured max deviation 2.2e-16; the test pins 1e-9).
Which candidate edges appear, with multi on Identical across two loads in one process, across 5 separate processes, and across 6 runs of krepp's own CLI
Which single edge is chosen, with multi=false No. Over 6 runs of krepp place --no-multi, 1 of 92 fragments alternated between two edges.

The multi=false case is worth understanding rather than just avoiding: krepp ranks candidates with a non-stable std::sort over a vector built by iterating a hash map keyed by tree-node pointer, then takes the last element. Equal-ranking candidates therefore land in an order that moves with the heap. The observed rate is low, but which read it hits is not something you can predict or pin.

If you need one deterministic row per read, leave multi on and pick with SQL (ORDER BY like_weight_ratio DESC, tie-broken on edge_num) rather than using multi=false. If you need reproducible like_weight_ratio values, round them — do not compare them for equality across runs.

Examples: 

-- Place reads from a FASTQ against a prebuilt index
CREATE TABLE reads AS
SELECT read_id, sequence1 FROM read_fastx('queries.fq');

SELECT * FROM place_krepp(query_table='reads', index_path='index_toy');

-- One row per read, chosen deterministically. Prefer this over multi=false,
-- which lets krepp break score ties on hash-map iteration order.
SELECT DISTINCT ON (fragment) fragment, edge_num, like_weight_ratio
FROM place_krepp(query_table='reads', index_path='index_toy')
ORDER BY fragment, like_weight_ratio DESC, edge_num;

-- An index built without a guide tree needs the backbone supplied
SELECT * FROM place_krepp(
    query_table='reads',
    index_path='index_notree',
    newick_path='backbone.nwk');

-- Feed placements into tree_resolve_placement, which names the same
-- quantities differently
CREATE TABLE placements AS
SELECT fragment AS fragment_id, edge_num AS edge_id,
       like_weight_ratio, distal_length, pendant_length
FROM place_krepp(query_table='reads', index_path='index_toy');

SELECT * FROM tree_resolve_placement('ref_tree', 'placements');

Availability: built when MIINT_ENABLE_KREPP is on, which excludes WASM and Windows. Check with SELECT * FROM miint_versions() WHERE library = 'krepp';, which reports the git describe --tags --dirty --always output that the build records for the submodule: v0.9.1-11-g3120503 for the current pin when the checkout has the v0.9.1 tag and the commits from it to the pin, 3120503 when it lacks either, and -dirty appended when the submodule has local changes. krepp's own version string for this pin is v0.10.0; upstream has no tag after v0.9.1.

Build a krepp index

Build the index place_krepp reads, from relations already in the database. The same linked krepp writes it and reads it, so the writer's version and the reader's version cannot drift apart — which matters more here than it sounds: an index written by a different krepp is misread rather than rejected (v0.9.0 swapped the recursion branches in Node::generate_tree), and as of 2026-09-10 bioconda's newest krepp is 0.9.1, older than the untagged v0.10.0 (31205033) this links.

Function signature:

krepp_index_create(sequence_table, output_path [, tree_table | newick_path, k, w, h, m, r, frac, sdust_t, sdust_w, threads])

Parameters: - sequence_table (VARCHAR, required, positional): Name of a table or view holding the reference sequences. Requires read_id and sequence1 columns, the same contract as the aligners and save_bowtie2_index; read_id may be VARCHAR, BIGINT or UUID, as it may for place_krepp's query_table.

Every distinct read_id becomes one krepp reference, and every row carrying that read_id becomes one record inside it. A multi-contig genome is therefore just several rows sharing a read_id — there is no separate mapping table. Reassemble chunked storage first, e.g. SELECT feature_id AS read_id, string_agg(chunk_data ORDER BY chunk_index) AS sequence1 FROM chunks GROUP BY feature_id. - output_path (VARCHAR, required, positional): Directory to write the index into. It may be absent, empty, or already hold other residues of the same index built with frac := false, which is how one build is split across jobs (see m, r and frac below). A different index, the same residue again, a second frac := true partial, an incomplete partial, or files with no complete partial beside them are refused, because krepp writes files whose names encode the resolved config and never clears what is already there. Files that are not index files — anything with an extension, or a name that does not start with cmer-, crecord-, inc-, metadata-, tree- or reflist- — are not inspected, so beside a complete partial they stay where they are and the build goes ahead. Nothing is deleted. - tree_table (VARCHAR): Name of a table or view holding the backbone tree in read_newick's row shape (node_index, parent_index, name, branch_length, edge_id). This is the form to prefer. The tree never becomes a file the caller has to render, so none of the strict-Newick hazards listed under place_krepp's newick_path apply — there is no indentation, no CRLF, no [&R] marker to get wrong. An edge_id column is carried through the build and comes back verbatim as place_krepp's edge_num. - newick_path (VARCHAR): Backbone tree as a Newick file, for callers who already have one. Held to exactly the same strictness place_krepp applies to its own newick_path, and checked before krepp sees it.

Exactly one of tree_table and newick_path is required. krepp will build without a backbone — it invents a tree from the reference names and then writes no tree file at all — but place_krepp cannot use that index without being handed a separate Newick, so the no-tree case is refused here rather than left as a trap. - k (INTEGER, default 29): k-mer length. krepp requires 19 ≤ k ≤ 31. - w (INTEGER, default k + 6): Minimizer window. Must be at least k. - h (INTEGER, default k − 16): Number of LSH positions. krepp requires 9 ≤ h ≤ 15, and h ≥ k − 16. - m (INTEGER, default 4), r (INTEGER, default 1), frac (BOOLEAN, default true): Partitioning of the k-mer space, exactly as krepp's --m, --r and --frac. These decide the index's filename suffix (-m4r1-frac).

r means different things depending on frac, and only one of them partitions. With frac := false, a k-mer is indexed when LSH(x) mod m == r, so each r is a disjoint slice. With frac := true — the default, matching krepp's own — the test is LSH(x) mod m <= r (ext/krepp/src/rqseq.cpp:133), making r a cumulative threshold: r := 2 contains everything r := 1 contains and more, and the partial is registered under every residue from 0 to r (ext/krepp/src/index.cpp:146-148). Measured on the 25-genome corpus in ext/krepp/test at m := 4 and the default k, total_num_kmers runs 8,675,085 / 13,491,823 / 17,366,632 for r = 1 / 2 / 3 — r := 3 already spans the whole space. krepp refuses to load two frac := true partials that differ in r (LSHF::check_compatible, ext/krepp/src/lshf.cpp:163-170), so a frac := true index is a single partial: krepp_index_create refuses to add a second one to an output_path, and place_krepp refuses a directory holding two.

Use frac := true to downsample a single index (m := 32, r := 16 keeps roughly half the k-mers), and frac := false to split one build across jobs. Only the second indexes a corpus too large for one pass. - sdust_t (INTEGER, default 0), sdust_w (INTEGER, default 0): SDUST low-complexity masking threshold and window. Both 0 disables masking; krepp warns that enabling it makes its subsampling model slightly inaccurate. - threads (INTEGER, default 1, maximum 256): Threads krepp spreads the build over. Available only when the extension was built with an OpenMP runtime present — SELECT * FROM miint_versions() WHERE library = 'krepp-openmp' returns a row when it was. On a build without one, threads := 1 is accepted and anything larger is refused rather than clamped, because a thread count that is quietly ignored is indistinguishable from a slow machine. The 256 ceiling is a guard rail, not a tuning limit: krepp passes the value to omp_set_num_threads unclamped, and when the runtime cannot create that many threads it aborts the process rather than failing the call.

Which files a threaded build changes: krepp merges each node's children in whatever order their tasks finish, so the colour arrays (cmer's subset column, and crecord) differ between a serial and a threaded build. The k-mers indexed and the bucket each lands in do not — inc-m4r1-frac was byte-identical across 1 and 4 threads, and total_num_kmers was unchanged. Note that no two builds are byte-identical across every file regardless of thread count, because metadata-*.txt records a build timestamp. Placement output is not a clean comparison in either direction — place_krepp is not reproducible run to run even single-threaded, for reasons that predate this parameter and are unrelated to it; the variation observed was confined to like_weight_ratio, with edge_num and the distances stable.

Scaling is sublinear and flattens early. Measured through krepp_index_create on a ten-core Apple M1 Max over the 25-genome corpus in ext/krepp/test, one run per point: about 10.1 s at one thread, 5.0 s at four, 4.6 s at eight.

Output schema: one row. - output_path (VARCHAR): The directory written. - k, w, h (INTEGER): The resolved parameters, read back from the index krepp wrote rather than echoed from the request — so a w or h left to default reports what was actually derived. - num_references (BIGINT): Distinct read_ids indexed. - status (VARCHAR): 'ok'. Failures raise.

Name matching is the failure this function exists to catch. krepp builds by walking the backbone tree and looking each tip up in its reference map, so the two name sets decide what ends up in the index:

Case krepp alone krepp_index_create
A reference the tree does not name Never visited. Absent from the index, with no message of any kind Raises, naming up to five of them
A tree tip with no reference Prints Genome skipped: to stderr, which nothing in a SQL session sees Warns with the count (see miint_warnings) and continues — a backbone wider than the reference set is legitimate

A reference name that Newick would have to quote (whitespace, or one of ()[]{},:;'") is also rejected: it would reach krepp quoted in the tree and bare in the map, and krepp matches the two by string equality, so the reference would silently contribute nothing. So is a name starting with > or @: krepp decides whether its input is a reference map or a FASTA from the first byte of the file, and the map's first line starts with a reference name.

Building one index across several jobs. Point every job at the same output_path, hold k, w, h, m and frac fixed, pass frac := false, and give each job a different r. krepp_index_create will add a partial to a directory that already holds others of the same index, and refuses to overwrite an r that is already there, to mix two hash configurations, or to add a partial whose k, w or h disagrees with what is on disk. place_krepp then reads the whole directory as one index. Jobs may run concurrently, but no two may be given the same r — the acceptance check is not atomic with krepp's write, so a duplicated residue is not caught.

Example:

-- References as rows: one row per contig, the genome id as read_id.
CREATE TABLE refs AS
SELECT regexp_extract(filepath, '([^/]+)\.fna$', 1) AS read_id, sequence1
FROM read_fastx('genomes/*.fna', include_filepath := true);

-- The backbone, with edge ids of our own that must survive to placement.
CREATE TABLE backbone AS SELECT * FROM read_newick('backbone.nwk');

SELECT * FROM krepp_index_create('refs', 'my_index', tree_table := 'backbone',
                                 k := 27, w := 35, h := 11);
-- my_index  27  35  11  25  ok

SELECT * FROM place_krepp(query_table := 'reads', index_path := 'my_index');

Behavior: - The build is single-threaded unless you ask otherwise. krepp parallelises its index build with OpenMP; whether those regions are compiled in depends on an OpenMP runtime being present at build time. Where one was, threads (above) spreads the build; where one was not, the pragmas are off and threads := 1 is the only accepted value. SELECT * FROM miint_versions() WHERE library = 'krepp-openmp' returns a row on a build that has them. - krepp writes progress to stderr — one line per tree node — which appears in the terminal rather than in the query result. - The index depends on the order rows arrive in. krepp's k-mer extraction for one reference is sensitive to the order of the records inside that reference's FASTA, and those records are written in the order sequence_table yields them — which SQL does not guarantee without an ORDER BY. Measured on a two-reference fixture, the same three rows grouped by read_id produce a k-mer table 8 bytes (one k-mer) larger than the same rows interleaved. The indexes are equally valid; they are not byte-identical. Add an ORDER BY to sequence_table if you need a reproducible artifact. - The function materialises one FASTA per reference in a temporary directory (under $TMPDIR) and hands krepp a name -> path map, because that is the only shape krepp's builder accepts with a guide tree: handed a single FASTA it switches to per-sequence mode, where each record is its own reference and a guide tree is refused outright. The temporary directory is removed when the call finishes, successfully or not. One file is open at a time, so the descriptor cost does not grow with the reference count; rows sharing a read_id that are not adjacent cost an extra open apiece, so grouping them pays. - Reference sequences are held to the IUPAC nucleotide alphabet (ACGTURYKMSWBDHVN, either case) — the same rule place_krepp applies to queries, for two reasons rather than one. A byte above 127 indexes past the end of krepp's 128-entry nucleotide table. And a newline followed by > inside a sequence would open a second FASTA record in that reference's file, which krepp neither rejects nor notices: it folds every record in a reference's file into the same tree leaf, so the reference would silently carry k-mers from content no row claimed. Offending rows raise, naming the byte and its offset. - A non-NULL sequence2 raises. A krepp reference is a single sequence; handed a paired-end relation, indexing only sequence1 would build a plausible index from half the data and report status = 'ok'. The column may exist as long as every value is NULL — the same contract the aligners' subject tables use. - Fatal errors inside krepp — an out-of-range k, a malformed tree, an unreadable input — are reported as SQL errors carrying krepp's own message. Without that, krepp's error_exit is a std::exit that would take the DuckDB process down with no error at all.

Availability: built when MIINT_ENABLE_KREPP is on, which excludes WASM and Windows. Check with SELECT 1 FROM duckdb_functions() WHERE function_name = 'krepp_index_create';.

Resolve placements

Resolve phylogenetic placements into a reference tree, returning a fully resolved tree with placed fragments as new tips. This exposes the insert_fully_resolved algorithm as a SQL-accessible table function.

Function signature:

tree_resolve_placement(tree_table, placements_table)

Parameters: - tree_table (VARCHAR): Name of a table or view containing tree data in read_newick schema (requires node_index and parent_index columns; name, branch_length, edge_id are optional) - placements_table (VARCHAR): Name of a table or view containing placement data: fragment_id (VARCHAR, BIGINT or UUID, as a read_id, so place_krepp's fragment of any id type can be used; the placed tip is named with its text), edge_id (BIGINT, INTEGER, UBIGINT or UINTEGER), and like_weight_ratio, distal_length and pendant_length (each DOUBLE or FLOAT)

Output schema: Same as read_newick (without filepath): - node_index (BIGINT): 0-based index of node in resolved tree - name (VARCHAR): Node label (placed fragments use their fragment_id as name) - branch_length (DOUBLE, nullable): Branch length - edge_id (BIGINT, nullable): Edge identifier (NULL for newly created nodes) - parent_index (BIGINT, nullable): Parent node's node_index (NULL for root) - is_tip (BOOLEAN): Whether node is a tip/leaf

Behavior: - Reads tree data from the tree table and builds a NewickTree - Reads placements and inserts each fragment as a new tip on the specified edge - Each placement creates 2 new nodes: an internal node (splitting the edge) and a fragment tip - Deduplicates placements by fragment_id, preferring the higher like_weight_ratio and, between ratios within 1e-9 of each other, the lower pendant_length; the kept placement is not always the one with the highest ratio (see Reproducibility below) - Multiple placements on the same edge are sorted by distal_length and inserted as a chain - Preserves original tip-to-tip distances in the tree - Works with both tables and views for either parameter - Schema is UNION ALL-compatible with read_newick

Reproducibility: The resolved tree is a function of the placement rows, not of the order they arrive in or of how they hash. Two things secure that.

The placements table is read in a pinned order (fragment_id, like_weight_ratio DESC, pendant_length, edge_id), because deduplication compares like_weight_ratio with a 1e-9 tolerance and "within tolerance" is not transitive — a ~ b and b ~ c do not imply a ~ c — so a single greedy pass over unordered rows can keep a different survivor each time. Note that this pins which row survives without making it the highest like_weight_ratio: a chain of rows each within tolerance of the next walks the pass down to the lowest pendant_length in the chain, which may sit outside tolerance of the maximum.

Edge grouping and chain order are then keyed on the data itself — edges in edge_id order, placements within an edge by distal_length then fragment_id's text — rather than on std::unordered_map iteration, whose order is unspecified and varies with insertion order and standard-library implementation. The fragment_id tiebreak matters because a placer reporting one position per edge leaves distal_length tied across every fragment sharing that edge; measured on place_krepp output, 735 of 735 edges had a single distinct distal_length.

Examples: 

-- Basic workflow: load tree, load placements, resolve
CREATE TABLE ref_tree AS
SELECT * FROM read_newick('reference.nwk');

CREATE TABLE placements AS
SELECT * FROM (VALUES
    ('seq1', 0::BIGINT, 0.95::DOUBLE, 0.05::DOUBLE, 0.001::DOUBLE),
    ('seq2', 1::BIGINT, 0.80::DOUBLE, 0.10::DOUBLE, 0.002::DOUBLE)
) AS t(fragment_id, edge_id, like_weight_ratio, distal_length, pendant_length);

SELECT * FROM tree_resolve_placement('ref_tree', 'placements');

-- Full jplace workflow: extract tree and placements from jplace file
CREATE TABLE jplace_tree AS
SELECT * FROM read_jplace_newick('results.jplace');

CREATE TABLE jplace_placements AS
SELECT fragment AS fragment_id, edge_num::BIGINT AS edge_id,
       like_weight_ratio, distal_length, pendant_length
FROM read_jplace('results.jplace');

-- Resolve and inspect the result
SELECT name FROM tree_resolve_placement('jplace_tree', 'jplace_placements')
WHERE is_tip = true
ORDER BY name;

-- Write resolved tree to Newick format
COPY (
    SELECT node_index, name, branch_length, edge_id, parent_index
    FROM tree_resolve_placement('ref_tree', 'placements')
) TO 'resolved.nwk' (FORMAT NEWICK);

-- Use a view for filtered placements
CREATE VIEW confident_placements AS
SELECT * FROM jplace_placements WHERE like_weight_ratio > 0.8;

SELECT COUNT(*) FROM tree_resolve_placement('jplace_tree', 'confident_placements')
WHERE is_tip = true;

Error conditions: - Tree table or placements table does not exist - Tree table missing required node_index or parent_index columns - Placements table missing required columns (fragment_id, edge_id, like_weight_ratio, distal_length, pendant_length), or a column of another type than listed above, when the query is planned and again when the placements are read - Placement references an edge_id not present in the tree - distal_length exceeds the edge's branch length - Negative distal_length or pendant_length

FastTree

Build an approximately-maximum-likelihood phylogenetic tree from a multiple sequence alignment using FastTree 2 (Price, Dehal, and Arkin 2010, Mol. Biol. Evol. 27:1641-1650). Returns one row per node (tips and internal nodes) with parent links suitable for joining back to the input or rendering to Newick via COPY ... (FORMAT 'newick').

Function signature: phylogeny_fasttree(table_name, [options])

Architecture: FastTree is statically linked into a separate process (gpl-boundary) that miint launches per query. The daemon is GPL-licensed; miint stays BSD-clean by communicating with it over POSIX shared memory and Arrow IPC. See docs/internals/embedded-tools.md for the protocol details and process-lifecycle invariants.

Requirements: - The gpl-boundary binary must be on PATH. Use phylogeny_fasttree_available() to check at runtime.

Parameters:

  • table_name (VARCHAR): Name of a table or view containing the input alignment. Must have columns name (VARCHAR) and sequence (VARCHAR), one row per tip. Minimum 3 sequences (FastTree cannot build a meaningful tree below that). All sequences should be the same length (a true MSA — pad with - if needed; see align_mafft).

The remaining parameters mirror the FastTree CLI flags. Defaults are the daemon's FastTree defaults; omitted parameters are not emitted into the wire config.

Sequence type and model:

  • seq_type (VARCHAR, default 'auto'): One of 'auto', 'nucleotide', 'protein'. Drives which substitution models are valid.
  • model (VARCHAR, default 'auto'): Substitution model. One of 'auto', 'jtt', 'lg', 'wag' (protein), 'jc', 'gtr' (nucleotide).
  • gtrrates (DOUBLE[], optional): GTR rate parameters (6 entries: AC, AG, AT, CG, CT, GT, all ≥ 0). Only valid with model='gtr'.
  • gtrfreq (DOUBLE[], optional): GTR base frequencies (4 entries: A, C, G, T, all ≥ 0). Only valid with model='gtr'.
  • gamma (BOOLEAN, optional): Use a discretized gamma distribution for among-site rate variation.
  • cat (BIGINT, optional, ≥ 1): Number of rate categories for CAT approximation.

Tree search:

  • nni (BIGINT, optional, ≥ 0): Number of NNI rounds. 0 disables NNI.
  • spr (BIGINT, optional, ≥ 0): Number of SPR rounds (subtree prune-regraft). 0 disables SPR.
  • mlnni (BIGINT, optional, ≥ 0): Number of ML NNI rounds.
  • mlacc (BIGINT, optional, ≥ 1): ML accuracy iterations per branch.
  • noml (BOOLEAN, optional): Skip the ML refinement phase entirely. Mutually exclusive with mlnni > 0.
  • slow (BOOLEAN, optional): Use slower / more accurate neighbor-joining (disables top hits).
  • bionj (BOOLEAN, optional): Use BIONJ joins. Mutually exclusive with nj.
  • nj (BOOLEAN, optional): Use vanilla NJ joins. Mutually exclusive with bionj.

Top-hits heuristics:

  • top (BOOLEAN, optional): Force top-hits heuristics on.
  • notop (BOOLEAN, optional): Force top-hits heuristics off.
  • topm (DOUBLE, optional, > 0): Top-hits multiplier (typical: 1.0-2.0). Only meaningful when top hits are enabled (rejected if notop=true or slow=true).

Bootstrap / support:

  • bootstrap (BIGINT, optional, ≥ 0): Number of resamples for SH-test branch support. 0 disables. Mutually exclusive with nosupport=true.
  • nosupport (BOOLEAN, optional): Suppress branch support output. Mutually exclusive with bootstrap > 0.

Pseudocounts (low-confidence trees):

  • pseudo (BOOLEAN, optional): Enable pseudocounts.
  • pseudo_weight (DOUBLE, optional, ≥ 0): Pseudocount weight. Requires pseudo=true.

Misc:

  • seed (BIGINT, optional): RNG seed. Pin for reproducibility (combine with threads=1).
  • threads (BIGINT, optional, ≥ 1): Thread count. Output is non-deterministic across thread counts — FastTree's parallel NJ/NNI/SPR sections produce floating-point variation. Force threads=1 for reproducibility.
  • verbose (BOOLEAN, optional): Forward FastTree's diagnostic output to the daemon log.
  • quote (BOOLEAN, optional): Quote tip names in the rendered Newick (when going through COPY FORMAT 'newick').

Cross-parameter validation rules (raised at SQL bind time, before any process spawn):

Constraint Error
bootstrap > 0 + nosupport=true "bootstrap=N cannot be combined with nosupport=true"
pseudo_weight without pseudo=true "pseudo_weight requires pseudo=true"
mlnni > 0 + noml=true "mlnni=N cannot be combined with noml=true"
Protein model + seq_type='nucleotide' "model='X' is a protein model and cannot be combined with seq_type='nucleotide'"
Nucleotide model + seq_type='protein' "model='X' is a nucleotide model and cannot be combined with seq_type='protein'"
gtrrates / gtrfreq without model='gtr' "gtrrates/gtrfreq are only valid with model='gtr'"
bionj=true + nj=true "bionj=true and nj=true are mutually exclusive"
top=true + notop=true "top=true and notop=true are mutually exclusive"
slow=true + top=true "slow=true disables top hits; cannot also set top=true"
topm + (notop=true or slow=true) "topm is only meaningful when top hits are enabled"
fastest=true Currently rejected — the embedded library API is a strict subset of the CLI's -fastest. Will be supported once the upstream surface catches up.

Output schema:

Column Type Description
node_index BIGINT 0-based row index, also the node identifier within this tree
parent_index BIGINT (nullable) node_index of this node's parent; NULL for the root
edge_id BIGINT (nullable) Edge identifier (matches the {N} placement notation, when used with read_jplace)
branch_length DOUBLE (nullable) Length of the edge from this node to its parent (NULL on the root)
support DOUBLE (nullable) SH-like branch support in [0,1] when bootstrap > 0; NULL otherwise
is_tip BOOLEAN True if this is a leaf (tip) node
name VARCHAR (nullable) Tip name (matches the input name column) for tips; NULL or empty for internal nodes
n_children BIGINT Number of direct children (0 for tips, 2 for binary internals, 3 for the unrooted root)

The schema is UNION ALL-compatible with read_newick if you select node_index, parent_index, name, branch_length, edge_id, is_tip (i.e., drop support and n_children).

Behavior:

  • One gpl-boundary daemon per query call (no cross-query connection caching yet).
  • Sequences are read from table_name via a fresh DuckDB connection (avoids re-entrant context.Query() deadlock — see docs/internals/reading-tables-views.md).
  • The full alignment is materialized in shared memory and submitted as a single Arrow IPC batch (FastTree is not a streaming algorithm).
  • The daemon is sent SIGTERM in the function's destructor; SIGKILL after a 30 × 10ms grace period if it doesn't exit cleanly.
  • FastTree's text output uses %.9f for branch lengths; miint's intermediate Arrow representation is full IEEE-754 double precision. Round-trip parity to the CLI binary is verified in test/sql/phylogeny_fasttree_parity.test.

Examples:

-- Basic workflow: align with MAFFT, then build a tree
CREATE TABLE seqs AS SELECT read_id AS name, sequence1 AS sequence FROM read_fastx('16s.fasta');

-- Default tree (auto-detected seq_type, FastTree defaults for everything else)
SELECT * FROM phylogeny_fasttree('seqs');

-- Reproducible run: pin seed and force threads=1
SELECT * FROM phylogeny_fasttree('seqs', seed := 42, threads := 1);

-- Nucleotide GTR with branch support from 100 bootstraps
SELECT * FROM phylogeny_fasttree(
    'seqs',
    seq_type := 'nucleotide',
    model := 'gtr',
    bootstrap := 100,
    seed := 42
);

-- GTR with explicit rate parameters
SELECT * FROM phylogeny_fasttree(
    'seqs',
    seq_type := 'nucleotide',
    model := 'gtr',
    gtrrates := [1.0, 2.5, 1.0, 1.0, 2.5, 1.0],
    gtrfreq  := [0.25, 0.25, 0.25, 0.25]
);

-- Render the tree to a Newick file
COPY (
    SELECT node_index, parent_index, name, branch_length, edge_id
    FROM phylogeny_fasttree('seqs', seed := 42, threads := 1)
) TO 'tree.nwk' (FORMAT 'newick');

-- Tip-only summary
SELECT name, branch_length AS terminal_branch_length
FROM phylogeny_fasttree('seqs', seed := 42)
WHERE is_tip
ORDER BY branch_length DESC;

-- Pipeline: align → tree → resolve placements (round-trip via read_newick)
COPY (
    SELECT node_index, parent_index, name, branch_length, edge_id
    FROM phylogeny_fasttree('seqs', seed := 42, threads := 1)
) TO 'ref.nwk' (FORMAT 'newick');

CREATE TABLE ref_tree AS SELECT * FROM read_newick('ref.nwk');
SELECT * FROM tree_resolve_placement('ref_tree', 'placements');

Error conditions:

  • gpl-boundary not on PATH (use phylogeny_fasttree_available() to gate)
  • Input table does not exist or is missing name / sequence columns
  • Fewer than 3 input sequences
  • Any of the cross-parameter rules above
  • Daemon process exits non-zero or returns a non-OK protocol response (exit message includes the daemon's stderr)

Reproducibility: With threads=1 and a fixed seed, the tree is bit-deterministic across runs of the same gpl-boundary build. With threads>1, FastTree's parallel sections produce floating-point variation in branch lengths and (rarely) topology — pin threads=1 for deterministic output.

Independent contrasts (PIC)

Compute Felsenstein's (1985) phylogenetic independent contrasts for one or more numeric per-tip traits over a tree. Per-tip traits are not statistically independent — related tips share evolutionary history — so naively correlating two traits across tips is biased. Independent contrasts is the standard correction: it transforms the correlated tip values into n-1 contrasts that are independent and identically distributed under a Brownian-motion model, which can then be correlated or regressed (through the origin) in downstream SQL.

Reference: Felsenstein, J. (1985) "Phylogenies and the Comparative Method." The American Naturalist 125(1):1–15.

Function signature:

phylo_independent_contrasts(tree_table, traits_table)

Parameters: - tree_table (VARCHAR): Name of a table or view containing tree data in read_newick schema (requires node_index and parent_index; name and branch_length are needed for a meaningful contrast; edge_id optional). The tree must be strictly bifurcating with finite, non-negative branch lengths (see requirements below). - traits_table (VARCHAR): Name of a table or view in long form with columns name, trait, value: - name: the tip identifier, matched to the tree's tip labels. Type-flexible — VARCHAR, UUID, and BIGINT are all accepted and matched to tip labels by their canonical text form (the same posture as read_id elsewhere), so a numeric or UUID feature key works directly. - trait (VARCHAR): the trait name; one trait is computed per distinct value. - value (DOUBLE): the numeric trait value for that tip.

Every trait must cover exactly the tree's tip set — no missing tips, no extra names, no NULL values. To analyze a subset of tips, shear_tree the tree to those tips first: pruning re-sums branch lengths, which changes the contrasts, so it cannot be skipped by simply omitting tips from the traits.

Output schema: long form, one row per (internal node, trait):

Column Type Description
node_index BIGINT The input tree's node_index for the internal node at which the contrast is computed
trait VARCHAR The trait name
contrast DOUBLE The standardized independent contrast (X_i − X_j) / sqrt(v_i + v_j)
ancestral_estimate DOUBLE The reconstructed (inverse-variance-weighted) trait value at the node
contrast_variance DOUBLE v_i + v_j, the contrast's expected variance (a property of the tree, identical across traits)

For a rooted bifurcating tree with n tips there are n-1 internal nodes and thus n-1 contrasts per trait.

Behavior / semantics: - Computed by Felsenstein's post-order pruning. At each internal node with children i, j and variance-extended branch lengths v_i, v_j: contrast = (X_i − X_j)/sqrt(v_i+v_j); the ancestral estimate is the inverse-variance-weighted mean (X_i·v_j + X_j·v_i)/(v_i+v_j); and the node's branch length is extended by (v_i·v_j)/(v_i+v_j) for use higher in the tree. - The sign of a contrast is arbitrary — it depends on which child is i vs j (here, node_index order). This is why the standard two-trait analysis is a regression through the origin over the paired contrasts (or a correlation after fixing one axis's sign). Results are reproducible: child order follows node_index. - contrast_variance is identical across traits for a given node, which is what makes multi-trait computation cheap (the tree is validated and its variances precomputed once). - The two-trait correlation / regression is intentionally not built in — layer it in SQL over the per-node output (example below).

Requirements and error conditions: - Strictly bifurcating. Every internal node must have exactly two children. A polytomy (> 2 children) errors and points to tree_resolve_multifurcations; a unifurcation (1 child) errors and points to shear_tree. Common gotcha: an unrooted Newick tree has a trifurcating root and must be resolved (or rooted) first. - Finite, non-negative branch lengths on every non-root edge. PIC divides by v_i + v_j, so an unspecified (NULL/NaN) or infinite length errors, and a negative length errors. The root's own branch length is unused. Zero-length internal edges are allowed (e.g. the connectors from tree_resolve_multifurcations); only a contrast whose variance v_i + v_j is exactly zero (both children zero-length) errors. - Trait coverage. A tip with no trait value, a trait value for a name that is not a tip, a duplicate (name, trait) pair, or a NULL value all error. - Tip names must be unique (they key the traits); tree_table / traits_table must exist and have the required columns.

Examples: 

-- Tree (must be bifurcating) and per-tip traits in long form.
CREATE TABLE tree AS SELECT * FROM read_newick('tree.nwk');
CREATE TABLE traits AS SELECT * FROM (VALUES
    ('sp1', 'body_mass', 3.2), ('sp2', 'body_mass', 1.1), ('sp3', 'body_mass', 2.7),
    ('sp1', 'metabolic_rate', 0.9), ('sp2', 'metabolic_rate', 0.4), ('sp3', 'metabolic_rate', 0.7)
) AS v(name, trait, value);

-- Contrasts for every trait.
SELECT * FROM phylo_independent_contrasts('tree', 'traits') ORDER BY trait, node_index;

-- Two-trait association via regression THROUGH THE ORIGIN over paired contrasts
-- (through the origin because each contrast's sign is arbitrary).
WITH c AS (SELECT * FROM phylo_independent_contrasts('tree', 'traits'))
SELECT sum(x.contrast * y.contrast) / sum(x.contrast * x.contrast) AS pic_slope
FROM c AS x JOIN c AS y USING (node_index)
WHERE x.trait = 'body_mass' AND y.trait = 'metabolic_rate';

-- If the tree has polytomies, resolve first (see the arbitrariness caveat there).
CREATE TABLE binary_tree AS SELECT * FROM tree_resolve_multifurcations('tree');
SELECT * FROM phylo_independent_contrasts('binary_tree', 'traits');

Cleanroom note: this is a from-scratch implementation of the algorithm as published in Felsenstein (1985); it does not derive from any GPL comparative-methods package.

Ancestral states (Brownian motion)

Reconstruct continuous ancestral states under a Brownian-motion (BM) model for one or more numeric per-tip traits — the maximum-likelihood estimate of each internal node's trait value, with a variance and 95% confidence interval. This is the sibling of independent contrasts (PIC): PIC's per-node ancestral_estimate is BM's post-order down-pass, which is the true ancestral value only at the root; full reconstruction adds a pre-order up-pass so every internal node conditions on the whole tree, not just the tips below it.

References: Felsenstein, J. (1985) The American Naturalist 125(1):1–15; Schluter, D. et al. (1997) "Likelihood of Ancestor States in Adaptive Radiation." Evolution 51(6):1699–1711.

Function signature:

phylo_ancestral_states(tree_table, traits_table)

Parameters: - tree_table (VARCHAR): a table/view in read_newick schema. Unlike PIC, the tree may be multifurcating (polytomies are handled directly — no need to tree_resolve_multifurcations). - traits_table (VARCHAR): long form name, trait, value — identical to PIC. name is type-flexible (VARCHAR/UUID/BIGINT, matched by canonical text); value is DOUBLE. Every trait must cover exactly the tree's tip set.

Output schema: long form, one row per (internal node, trait):

Column Type Description
node_index BIGINT The input tree's node_index for the internal node
trait VARCHAR The trait name
estimate DOUBLE The marginal ML ancestral estimate (mean) at the node
variance DOUBLE Variance of the estimate (σ̂² × the node's structural variance)
ci_low DOUBLE 95% CI lower bound: estimate − 1.959964·sqrt(variance)
ci_high DOUBLE 95% CI upper bound: estimate + 1.959964·sqrt(variance)

Behavior / semantics: - Computed by Gaussian belief propagation on the tree: a post-order down-pass (each internal node is the precision-weighted mean of its children's messages — this reduces to PIC at a bifurcation but is correct for any arity) plus a pre-order up-pass (an O(n) leave-one-out that folds in the rest of the tree). The result is exact on a tree. - The root estimate equals PIC's root ancestral_estimate (a free cross-check); interior estimates generally differ from PIC's down-pass values because the up-pass adds information from outside each subtree. - The BM rate σ̂² is estimated by REML (sum of squared standardized contrasts divided by n_tips − 1, the unbiased/PIC-consistent estimator). The alternative full-ML divisor n_tips would differ by a factor (n_tips−1)/n_tips. variance and the CI are σ̂² times each node's rate-independent structural variance; z = 1.959964 = qnorm(0.975).

Requirements and error conditions: - Multifurcations are supported; a unifurcation (internal node with a single child) errors and points to shear_tree (it carries no information). - Finite, non-negative branch lengths on every non-root edge, and a strictly positive length on every tip edge. A zero-length tip edge would pin the estimate and break the leave-one-out, so it errors (a deliberate difference from PIC, which tolerates one zero-length tip per contrast). Zero-length internal edges are allowed. - Trait coverage (missing tip, extra name, duplicate (name, trait), NULL value) errors as in PIC. A non-finite (NaN/infinite) trait value also errors — it would silently poison the shared rate and the whole output, so it fails loud.

Examples: 

CREATE TABLE tree AS SELECT * FROM read_newick('tree.nwk');   -- may be multifurcating
CREATE TABLE traits AS SELECT * FROM (VALUES
    ('sp1', 'body_mass', 3.2), ('sp2', 'body_mass', 1.1), ('sp3', 'body_mass', 2.7)
) AS v(name, trait, value);

-- Ancestral estimate + 95% CI at every internal node.
SELECT * FROM phylo_ancestral_states('tree', 'traits') ORDER BY trait, node_index;

-- Cross-check: the root estimate equals PIC's root ancestral_estimate.
SELECT a.estimate, p.ancestral_estimate
FROM phylo_ancestral_states('tree', 'traits') a
JOIN phylo_independent_contrasts('tree', 'traits') p USING (node_index, trait)
WHERE a.node_index = (SELECT node_index FROM tree WHERE parent_index IS NULL);

Cleanroom note: a from-scratch implementation of BM ancestral reconstruction from the published algorithm (Felsenstein 1985; Schluter et al. 1997). It does not derive from any GPL comparative-methods package. Correctness is validated against an independent generalized-least-squares (phylogenetic VCV) ground truth, which is exact on a tree. (We don't use ape::ace(method="ML") as the reference: its point estimates agree with this GLS truth, but it reports the BM rate σ² on a different scale, so its variances/CIs aren't directly comparable.)

Ancestral parsimony (Sankoff)

Reconstruct discrete ancestral states by parsimony — the state assignments that minimize the total number (or cost) of state changes over the tree — for one or more categorical per-tip traits. The engine is Sankoff's dynamic program, which is correct for any tree shape and any substitution-cost matrix; with the default unit cost (0 for no change, 1 for any change) it is Fitch parsimony.

References: Sankoff, D. (1975) SIAM J. Appl. Math. 28(1):35–42; Fitch, W.M. (1971) Systematic Zoology 20(4):406–416.

Function signature:

phylo_ancestral_parsimony(tree_table, traits_table [, cost_matrix_table])

Parameters: - tree_table (VARCHAR): a table/view in read_newick schema. Branch lengths are ignored — topology-only trees are fine. Multifurcations and unifurcations are both supported. - traits_table (VARCHAR): long form name, trait, value. name is type-flexible as elsewhere; value is the categorical state, cast to VARCHAR (so integer- or text-coded states both work). Every trait must cover exactly the tree's tip set. - cost_matrix_table (VARCHAR, optional): long form from_state, to_state, cost (a DOUBLE ≥ 0, the cost of a change from from_state to to_state). When supplied it defines the state alphabet (the sorted union of from_state/to_state), which may include states not observed at any tip (a Sankoff intermediate). Every off-diagonal ordered pair over that alphabet must be present (a missing transition cost errors — there is no sensible default); the diagonal defaults to 0. Omit it for the unit-cost (Fitch) default, whose alphabet is that trait's sorted distinct observed states.

Output schema: long form, one row per (internal node, trait, state):

Column Type Description
node_index BIGINT The input tree's node_index for the internal node
trait VARCHAR The trait name
state VARCHAR A state from the trait's alphabet
in_mpr BOOLEAN Whether this state is in the node's most-parsimonious-reconstruction (MPR) set
min_cost DOUBLE The minimum total tree cost with this node fixed to this state

Behavior / semantics: - A post-order Sankoff down-pass computes, per node/state, the min cost of the subtree below it; a pre-order pass computes the rest-of-tree cost; their sum is min_cost (the min total cost with the node pinned to that state). The whole-tree parsimony score is min(min_cost), and in_mpr is true exactly when min_cost equals that score. - Ties are first-class: an ambiguous node simply has more than one in_mpr = true state. Filter with WHERE in_mpr for the MPR set; the non-MPR rows give the extra cost of each alternative state. - Multifurcation and unifurcation are both handled by the same recurrence (arity-agnostic).

Requirements and error conditions: - Tip names must be unique; trait coverage (missing tip, extra name, duplicate (name, trait), NULL) errors. - A cost_matrix_table with a missing off-diagonal entry, a non-zero diagonal, a duplicate pair, a negative/non-finite cost, or a tip whose observed state is absent from the matrix alphabet all error.

Examples: 

CREATE TABLE tree AS SELECT * FROM read_newick('tree.nwk');
CREATE TABLE traits AS SELECT * FROM (VALUES
    ('sp1', 'habitat', 'marine'), ('sp2', 'habitat', 'marine'), ('sp3', 'habitat', 'terrestrial')
) AS v(name, trait, value);

-- Fitch (unit-cost) reconstruction; the MPR state(s) per node:
SELECT node_index, state, min_cost
FROM phylo_ancestral_parsimony('tree', 'traits') WHERE in_mpr ORDER BY node_index;

-- Weighted Sankoff: an ordered character 0 - 1 - 2 where a direct 0<->2 change is dear.
CREATE TABLE cost AS SELECT * FROM (VALUES
    ('0','1',1.0), ('1','0',1.0), ('1','2',1.0), ('2','1',1.0), ('0','2',3.0), ('2','0',3.0)
) AS v(from_state, to_state, cost);
SELECT * FROM phylo_ancestral_parsimony('tree', 'traits', 'cost') ORDER BY node_index, state;

Cleanroom note: a from-scratch implementation of Sankoff (1975) parsimony (Fitch 1971 for the unit-cost case). It does not derive from any GPL comparative-methods package. Correctness (including multifurcations and arbitrary cost matrices) is validated against an independent brute-force enumerator over all state assignments.

Ancestral maximum likelihood (Mk)

Reconstruct discrete ancestral states by maximum likelihood under the Mk model (a k-state continuous-time Markov model), returning the marginal posterior probability distribution over states at each internal node. Unlike parsimony, this uses branch lengths and reports a full probability distribution and a model fit.

Three models are supported: ER (equal rates — all state changes share one rate), SYM (symmetric rates — q_ij = q_ji, with k(k−1)/2 independent off-diagonal rates), and ARD (all rates different — every one of the k(k−1) off-diagonal rates is free). ER and SYM are reversible; ARD is not (its reconstruction depends on where the tree is rooted). All three use a uniform root prior 1/k.

References: Lewis, P.O. (2001) Systematic Biology 50(6):913–925 (Mk model); Felsenstein, J. (1981) J. Mol. Evol. 17(6):368–376 (pruning likelihood); Yang, Z., Kumar, S. & Nei, M. (1995) Genetics 141(4):1641–1650 (marginal reconstruction); Nelder, J.A. & Mead, R. (1965) Computer Journal 7(4):308–313 (simplex optimizer, for the SYM/ARD fits); Higham, N.J. (2005) SIAM J. Matrix Anal. Appl. 26(4):1179–1193 (scaling-and-squaring matrix exponential, for the ARD P(t)).

Function signature:

phylo_ancestral_ml(tree_table, traits_table, model := 'ER', rate := NULL)

Parameters: - tree_table (VARCHAR): a table/view in read_newick schema. Branch lengths are used (the transition probabilities depend on them). Multifurcations and unifurcations are supported. - traits_table (VARCHAR): long form name, trait, value (state cast to VARCHAR), as for parsimony. The per-trait alphabet is its sorted distinct observed states. - model (VARCHAR, named, default 'ER'): the substitution model — 'ER' (equal rates), 'SYM' (symmetric rates), or 'ARD' (all rates different). Any other value errors. - rate (DOUBLE, named, optional): applies to ER only. If given (a finite positive number), the single rate is fixed to it; otherwise it is fitted by maximum likelihood. It has no meaning for SYM/ARD (which fit a whole rate matrix) — supplying rate with model := 'SYM' or 'ARD' errors.

Output schema: long form, one row per (internal node, trait, state):

Column Type Description
node_index BIGINT The input tree's node_index for the internal node
trait VARCHAR The trait name
state VARCHAR A state from the trait's alphabet
probability DOUBLE Marginal posterior probability that the node is in this state given the data; the per-node probabilities sum to 1
rate DOUBLE ER: the fitted (or fixed) scalar rate, the same value on every row of the trait. SYM/ARD: NULL (a rate matrix has no single scalar rate)
log_likelihood DOUBLE The model log-likelihood at the fitted/fixed parameters — the same value on every row of the trait

Behavior / semantics: - Computed by Felsenstein's pruning algorithm with per-node rescaling (numerically stable on deep trees), a uniform root prior (the stationary distribution 1/k, shared by ER and SYM), and a marginal two-pass so every internal node conditions on the whole tree. - ER: the transition matrix is closed-form (P(t)[i][j] = 1/k + (δ_ij − 1/k)·e^{−k·rate·t}), so no eigendecomposition is needed and it scales to large trees. The single rate is fitted by a grid-scan + golden-section search on the log-likelihood. - SYM: the transition matrix P(t) = exp(Q·t) is built from a symmetric eigendecomposition of Q (P(t) = U·diag(e^{λt})·Uᵀ); the k(k−1)/2 off-diagonal rates are fitted jointly by a Nelder–Mead simplex search (warm-started from the ER fit, so SYM's log-likelihood is always at least ER's). At k = 2 SYM has one free rate and reduces exactly to ER. - ARD: Q is a general (non-symmetric) generator, so P(t) = exp(Q·t) is built by a general matrix exponential (scaling-and-squaring). The k(k−1) rates are fitted by a multi-start Nelder–Mead search (the ER warm start plus several deterministic fixed-seed restarts, so the fit is reproducible). This is a best-effort optimizer: on ARD's high-dimensional, multimodal surface it may settle at a local optimum, and gradient-based tools can reach boundary (rate→0) solutions this log-space search does not. Because ARD is non-reversible, the reconstruction depends on the root (there is no pulley-principle invariance). ARD does not scale to large state counts (k(k−1) simplex parameters) or, as fast as ER, to very large trees. - State-count caps: SYM rejects k > 8 and ARD rejects k > 6 (the simplex fit becomes unreliable) — use ER for high-state traits. - Per-trait rate (ER only; NULL for SYM/ARD) and log_likelihood are surfaced (denormalized onto every row) so a SELECT DISTINCT trait, rate, log_likelihood recovers the fit. - ER and SYM output match ape::ace(type="discrete", model=...) posteriors closely — to ~1e-8 (ER) and ~1e-6 (SYM), i.e. optimizer/iteration tolerance — and the ER rate to ~1e-8. For ARD the reported fit is not guaranteed to be the global optimum: the search is in log-rate space, so it cannot reach boundary solutions where a rate is exactly 0 (a common true MLE for sparse discrete data with rare or absent transitions) — a box- or gradient-constrained optimizer (e.g. ape) may find a higher-likelihood fit on the same data. ARD's machinery (pruning, P(t), the marginal up-pass) is nonetheless validated exactly against an independent brute-force likelihood at whatever rates are fitted; only the optimizer, not the reconstruction math, is best-effort.

Requirements and error conditions: - Finite, non-negative branch lengths on every non-root edge, and a strictly positive length on every tip edge (as in ancestral states); zero-length internal edges are allowed. A branch so small that e^{−k·rate·t} rounds to exactly 1 (an effectively-zero edge) can make the likelihood degenerate and errors (rather than emitting a NaN/−inf silently). - Tip names unique; trait coverage errors as elsewhere. An unsupported model, a non-positive / non-finite rate, a rate combined with model := 'SYM'/'ARD', or a state count above a model's cap (SYM k > 8, ARD k > 6), errors.

Examples: 

CREATE TABLE tree AS SELECT * FROM read_newick('tree.nwk');
CREATE TABLE traits AS SELECT * FROM (VALUES
    ('sp1', 'diet', 'carnivore'), ('sp2', 'diet', 'herbivore'), ('sp3', 'diet', 'herbivore')
) AS v(name, trait, value);

-- Fit the ER rate by ML; posterior distribution at each internal node.
SELECT node_index, state, probability FROM phylo_ancestral_ml('tree', 'traits') ORDER BY node_index, state;

-- The fitted rate and model fit per trait.
SELECT DISTINCT trait, rate, log_likelihood FROM phylo_ancestral_ml('tree', 'traits');

-- Reconstruct at a fixed (known) rate instead of fitting (ER only).
SELECT * FROM phylo_ancestral_ml('tree', 'traits', rate := 0.5);

-- Symmetric-rates (SYM) or all-rates-different (ARD): the full rate matrix is fitted; the
-- `rate` column is NULL. ARD additionally distinguishes gain vs loss (direction) of a state.
SELECT node_index, state, probability FROM phylo_ancestral_ml('tree', 'traits', model := 'SYM') ORDER BY node_index, state;
SELECT node_index, state, probability FROM phylo_ancestral_ml('tree', 'traits', model := 'ARD') ORDER BY node_index, state;

Cleanroom note: a from-scratch implementation of Mk maximum-likelihood reconstruction from the published algorithms (Lewis 2001; Felsenstein 1981; Yang, Kumar & Nei 1995; Nelder & Mead 1965 for the SYM/ARD optimizer; Higham 2005 for the ARD matrix exponential). It does not derive from any GPL comparative-methods package. The P(t) transition matrices use Eigen (MPL-2.0, a permissive license compatible with this project's Modified-BSD license; EIGEN_MPL2_ONLY is defined): a SelfAdjointEigenSolver for SYM and the MatrixExponential (scaling-and-squaring Padé) for ARD. Correctness is validated against an independent brute-force likelihood over all state assignments (with an independent matrix exponential for SYM/ARD). The ER and SYM conventions are confirmed to match ape::ace(type="discrete", model=...) numerically (that tool is used offline as an oracle only — its code is never read, linked, or distributed).