SplitAligner, end to end

Most comparative-genomics questions are answered by comparing a number attached to a branch — a rate, a residual, a path length — across many genes. Did this lineage speed up? Did constraint relax here? The comparison only means something if that branch denotes the same evolutionary object in every gene that contributes to it.

In practice it often does not. Genes differ in which taxa they recover, and gene trees do not always agree with the species tree. When a taxon is missing, two adjacent branches of the reference tree can collapse into one. The merged branch has a perfectly correct length — the sum of its parts — but it is no longer either of the branches that formed it.

A pipeline that files that number under one of them produces a value that is numerically right and biologically misassigned. Nothing crashes, no NA appears, and no check that only verifies arithmetic will catch it. The practical consequence is that a gene can appear to evolve unusually fast or slow on a lineage for reasons that have nothing to do with evolution on that lineage — and without a record of which cells are affected, there is no way to tell those cases from real ones.

SplitAligner's response is to make branch identity explicit, and to give every way a branch can fail to be observable its own name. The animated explainer makes the same argument visually; the rest of this page is how to run it.

SplitAligner identifies a branch not by its drawn position or by node order, but by the bipartition (split) it induces on the taxa. For a gene g that samples only some taxa T_g, each species-tree branch b is projected onto the taxa actually present:

On the topology-congruent species-tree restriction, two non-deleted primitive branches map to the same reduced-tree edge exactly when their projected splits coincide. They do not become the same original branch: they become indistinguishable for that gene's taxon set, which is a fusion. If a projected side is empty, the branch instead has a structural absence. For an empirical free-topology gene tree, SplitAligner separately asks whether each independently observable projected split is recovered.

SplitAligner's algorithm is representation-independent. Because identity lives in the split, how a tree is written does not matter: every edge of a gene tree is reduced to a canonical unrooted split before any comparison. A gene tree may therefore be rooted or unrooted, binary or multifurcating, and written in any child order; rooted input is simply treated as unrooted: the two half-edges at a degree-2 root become one split; their lengths are summed when both are finite numeric values, and otherwise the split carries no numeric length. Its topology may be constrained to the species tree (fix) or completely free. Two practical requirements remain: the taxa of each gene tree must be a subset of the species-tree taxa, with names matching exactly, and the gene trees should carry branch lengths, since those are the values that get filed into the matrix. This is verified by a serialization- and input-variant regression battery on the bundled example: equivalent inputs give identical matrices, byte for byte when the species tree is unchanged and cell for cell after mapping the B aliases when it is re-serialized.

# Split-based branch-mapping algorithm

# Stage 1 — coverage level: uses only S and T_g, never the gene tree
for each gene g:
    T_g ← taxa of the gene tree (a subset of S)               # Step 1
    for each species-tree branch b:                           # Steps 2–3
        σ_g(b) ← canon(A_b ∩ T_g | B_b ∩ T_g)                 # projection
        if one side of σ_g(b) is empty:
            state(g,b) ← NA_struct                            # no projected identity
        else if b is internal and one side of σ_g(b) is a single taxon x:
            endpoint-collapsed: σ_g(b) coincides with the terminal branch of x,
            so b is kept for fusion bookkeeping, never as an independent branch
    group the retained branches by σ_g(b)                     # Step 4
        a group with ≥ 2 members is a fusion group: its members are
        indistinguishable on T_g and share one composite coordinate b|b'
    every retained (g,b) now has a settled identity; its value is still unknown

# Stage 2 — gene-tree level: each retained projected split consults the gene tree once
    Σ(G_g) ← canonical unrooted splits of G_g                 # Step 5
    for each retained projected split σ and its group:
        if σ ∉ Σ(G_g):            generic NA for every member   # cause not yet named
        else if the group has ≥ 2 members:
                                NA_fuse for every member      # the length is filed on b|b'
        else:                   Mapped — branch length recorded

# Stage 3 — finalize: paired fixed/free evidence names the remaining NAs
    NA_fuse stands only if the composite coordinate b|b' carries a numeric length
    a free-side generic NA becomes NA_topo only when the fixed-side primitive
    coordinate is Mapped with a numeric length; otherwise it stays residual NA

The workflow has two modes: matrix builds the coordinate system and projects every gene onto it; finalize then classifies the resulting missingness where the paired evidence permits.

You need Perl 5 (already present on macOS and most Linux systems). SplitAligner uses only core modules — no CPAN packages to install.

That is the entire installation. The main controller is SplitAligner.pl; the per-step helpers live in scripts/ and are called for you.

Species tree — one Newick tree. Branch lengths and support values are allowed and simply ignored during split-based mapping.

Gene trees — Newick, one record per line, each line starting with a gene identifier immediately followed by its tree. Taxon names must match the species tree.

SplitAligner distinguishes two kinds of gene-tree input, and you typically supply both:

• free — each gene tree inferred independently (its topology may disagree with the species tree);
• fix — each gene tree inferred under the species-tree topology (branch lengths only).

Comparing the two is what later separates coverage-driven absence from discordance-driven absence.

The repository ships a small example under examples/302mammal/ (a subset of free and fixed gene trees) for a fast end-to-end smoke test. From the repository root:

This runs the three program calls below in order. The rest of the tutorial unpacks what each one does and produces.

Run matrix once per gene-tree set. It labels the species-tree branches, projects every gene onto them, and writes a gene × branch matrix.

Key outputs (prefixed by --label):

• species_tree.splits.txt — the canonical branch coordinate system (B1, B2, … on the species-tree backbone). This is the shared axis every gene is measured against.
• <label>.matrix_no_fuse.txt — primitive branches only.
• <label>.matrix_with_fuse.txt — the same, plus explicit fused-branch columns such as B12|B47 for branches that merge under taxon pruning.
• <label>_splits/ and <label>_split_branch_label/ — per-gene split representations and their mapping onto the backbone.

A cell holds a numeric branch length where the branch's projected split was recovered, or NA where it was not — yet. Naming those NAs is the next step.

Give finalize both matrices. It first identifies primitive NAs represented by a numeric fused coordinate (NA_fuse), then compares the free and fixed matrices on shared genes to separate structural absence (NA_struct) from topology-induced absence (NA_topo). If fixed-side numeric evidence is unavailable, a residual generic NA is retained. Passing --species_tree also computes branch-wise support.

Outputs:

• final.free.na_classified.txt / final.fix.na_classified.txt — classified matrices, with residual generic NA retained when the fixed side lacks numeric evidence;
• <input>.na_fuse.txt — the intermediate fusion-rescue tables;
• final.support_b.txt and species_tree.support_b.nwk — branch-wise support (see step 7).

The comparison is defined only for genes present in both inputs; if there are no shared genes, SplitAligner stops with an error rather than guessing.

Every gene × branch cell ends up as exactly one of the values below. A residual generic NA is an intentional bookkeeping state: it is not promoted to NA_topo without fixed-side numeric evidence.

The categories partition the axis — they add up to the whole ledger:

Once NA is split into why, downstream rate or selection analyses can handle, filter, and audit the states explicitly instead of treating every absence as the same opaque value.

When you pass --species_tree, finalize writes final.support_b.txt: for each backbone branch, how much numeric evidence it retains in the free trees relative to the fixed trees, on shared genes.

The companion species_tree.support_b.nwk is a standard Newick tree with each internal branch's Support written in the bootstrap position, so you can open it in any tree viewer (e.g. FigTree) and see concordance mapped straight onto the topology.

The same two modes run at full scale. The bundled preprint example uses 2,275 free-topology and 2,275 fixed-topology gene trees across 302 mammals:

On the shared concordance subproblem, SplitAligner's Support(b) closely matches PhyParts (Pearson r = 0.99997 across 299 internal branches). On the reported single-machine benchmark, the full SplitAligner workflow took 9 min 52 s versus 1 h 50 min 53 s for PhyParts in its quick concordance setting — about 11× faster on this dataset. This is a calibration and runtime comparison, not a claim of functional equivalence.

Let taxon pruning act alone, with no gene-tree discordance. The full B1–B17 axis stays frozen while each deletion updates the graph-determined ledger: observed, NA_fuse, or NA_struct. In this experiment NA_topo is always zero.