SplitAligner：从头到尾走一遍

比较基因组学的大多数问题，都是靠比较「挂在某根分枝上的一个数值」——速率、残差、路径长度——在许多基因之间的差异来回答的。这个系统的速率上升了吗？这里的约束放松了吗？只有当那根分枝在每一个参与比较的基因里都指同一个进化对象时，这种比较才有意义。

实际情况往往不是这样。不同基因回收到的类群不同，基因树也未必与物种树一致。当某个类群缺失时，参考树上相邻的两根分枝可能坍缩成一根。融合后那根分枝的长度完全正确——它是各组成部分之和——但它已经不再是形成它的任何一根分枝了。

把这个数值记在其中一根名下的流程，产生的是一个数值正确、而生物学归属错误的值。程序不会崩溃，不会出现 NA，任何只检查算术是否正确的检验都抓不到它。实际后果是：一个基因在某个系统上看起来异常快或异常慢，原因可能与那个系统上的进化毫无关系；而如果没有一份「哪些单元格受了影响」的记录，就无法把这类情况和真实情况区分开。

SplitAligner 的应对是把分枝身份显式化，并为分枝失去可观测性的每一种方式各起一个名字。动画讲解用视觉方式讲同一个论证；本页余下部分讲的是怎么运行它。

SplitAligner 识别一根枝，靠的不是它画在哪里、也不是节点顺序，而是它在类群上诱导出的二分划分（split）。对于只采样了部分类群 T_g 的基因 g，把物种树的每根枝 b 投影到实际存在的类群上：

在与物种树限制拓扑一致的基准层中，两根未被删除的原始枝投影划分相同，当且仅当它们落到同一条约化枝上。这并不是说两根原始枝变成了同一根；而是它们在该基因的类群集合下无法再独立区分，形成融合。若投影后有一侧为空，则是结构性缺失。对于经验 free 拓扑基因树，SplitAligner 还会另行判断每个可独立观测的投影划分是否被找回。

SplitAligner 的算法是表示非依赖的。既然枝的身份就是它的划分，一棵树怎么写就无关紧要：比较之前，基因树的每一条枝都会先归约为规范化的无根划分。因此基因树可以有根或无根、二叉或多叉、子节点顺序任意；有根的输入一律按无根树处理：二度根处的两条半枝合并为同一个划分；两者都是有限数值时枝长相加，否则该划分不带数值枝长。拓扑可以约束在物种树上（fix），也可以完全自由。实际使用只有两点要求：每棵基因树的类群必须是物种树类群的子集，名称完全一致；基因树应带枝长，因为枝长正是要填进矩阵的值。这一点已由自带示例上的序列化与输入变体回归测试验证：等价输入给出相同的矩阵，物种树不变时逐字节相同，物种树换一种写法时经 B 别名映射后逐格相同。

# Split-based branch-mapping algorithm

# 第一层 —— 覆盖层：只用 S 和 T_g，完全不看基因树
for 每个基因 g:
    T_g ← 基因树的类群（S 的子集）                               # 步骤 1
    for 每根物种树枝 b:                                     # 步骤 2–3
        σ_g(b) ← canon(A_b ∩ T_g | B_b ∩ T_g)                 # 投影
        if σ_g(b) 有一侧为空:
            state(g,b) ← NA_struct                            # 没有投影身份
        else if b 是内部枝且 σ_g(b) 有一侧只剩一个类群 x:
            端点塌缩：σ_g(b) 与 x 的终端枝划分相同，
            所以 b 只做融合记账，永远不作为独立枝
    把保留下来的枝按 σ_g(b) 分组                                # 步骤 4
        成员 ≥ 2 的组是融合组：成员在 T_g 上无法区分，
        共用一个复合坐标 b|b'
    此时每个保留的 (g,b) 身份已定，值尚未知

# 第二层 —— 基因树层：每个保留的投影划分查一次基因树
    Σ(G_g) ← G_g 的规范无根划分                               # 步骤 5
    for 每个保留的投影划分 σ 及其分组:
        if σ ∉ Σ(G_g):            所有成员记 generic NA        # 原因尚未命名
        else if 组内成员 ≥ 2:
                                所有成员记 NA_fuse              # 枝长记在 b|b' 上
        else:                   Mapped —— 记录枝长

# 第三层 —— finalize：用配对的 fixed/free 证据给剩下的 NA 命名
    NA_fuse 只在复合坐标 b|b' 带有数值枝长时成立
    free 侧的 generic NA 只有在 fixed 侧原始坐标为 Mapped 且带数值时
    才成为 NA_topo；否则保持为 residual NA

整个流程有两个模式：matrix 构建坐标系并把所有基因投影上去；finalize 再根据成对证据区分缺失状态。

你需要 Perl 5（macOS 和大多数 Linux 自带）。SplitAligner 只用核心模块——无需安装任何 CPAN 包。

安装到此为止。主控脚本是 SplitAligner.pl；各步骤的辅助脚本在 scripts/ 中，会被自动调用。

物种树 —— 一棵 Newick 树。允许带枝长和支持度，在基于划分的映射中会被直接忽略。

基因树 —— Newick 格式，一行一条记录，每行以基因标识符开头，紧接着是该基因的树。类群名必须与物种树一致。

SplitAligner 区分两类基因树输入，通常两者都要提供：

• free —— 每个基因独立推断的树（拓扑可以与物种树不一致）；
• fix —— 在物种树拓扑约束下推断的树（只有枝长）。

比较这两者，正是之后区分覆盖度导致的缺失与不一致导致的缺失的关键。

仓库在 examples/302mammal/ 下附带了一个小示例（free / fix 基因树的子集），用于快速的端到端冒烟测试。在仓库根目录下运行：

它会依次运行下面三次程序调用。教程余下部分逐一拆解每一步做了什么、产出什么。

每组基因树运行一次 matrix。它给物种树的枝打标签，把所有基因投影上去，写出一张基因 × 枝矩阵。

主要输出（以 --label 为前缀）：

• species_tree.splits.txt —— 标准的枝坐标系（物种树骨架上的 B1、B2……）。这是所有基因据以度量的共同坐标轴。
• <label>.matrix_no_fuse.txt —— 仅原始枝。
• <label>.matrix_with_fuse.txt —— 在此基础上，额外加入显式的融合枝列（如 B12|B47），对应因类群修剪而合并的枝。
• <label>_splits/ 和 <label>_split_branch_label/ —— 每个基因的划分表示，及其到骨架的映射。

若枝的投影划分被找回，单元格即为其枝长（数值）；否则为 NA——但只是暂时的。为这些 NA 命名，就是下一步。

把两张矩阵都交给 finalize。它先识别由数值融合坐标解释的原始枝 NA，标为 NA_fuse；再在共享基因上比较 free 与 fix 矩阵，把结构性缺失（NA_struct）与拓扑导致的缺失（NA_topo）分开。若 fix 侧没有可供判断的数值证据，则保留残余 NA。传入 --species_tree 时，还会计算枝级支持度。

输出：

• final.free.na_classified.txt / final.fix.na_classified.txt —— 完成状态分类的矩阵；fix 侧缺少数值证据时会保留残余 NA；
• <input>.na_fuse.txt —— 融合救回的中间表；
• final.support_b.txt 和 species_tree.support_b.nwk —— 枝级支持度（见第 7 步）。

比较只对两个输入都存在的基因有定义；若没有共享基因，SplitAligner 会直接报错停止，而不去猜。

每个基因 × 枝单元格最终都恰好取下列一种值。残余 NA 是有意保留的记账状态：没有 fix 侧数值证据时，不能把它升级为 NA_topo。

这些类别对坐标轴构成一个划分——加起来就是整本台账：

一旦 NA 按原因拆开，下游的速率或选择分析就能明确地处理、过滤和审计各类状态，而不是把所有缺失都压成同一个不透明的值。

当你传入 --species_tree，finalize 会写出 final.support_b.txt：对骨架上的每根枝，给出它在共享基因上、free 树相对 fix 树保留了多少数值证据。

配套的 species_tree.support_b.nwk 是一棵标准 Newick 树，把每根内部枝的 Support 写在自展值（bootstrap）位置上，可在任意树浏览器（如 FigTree）中打开，直接把一致度映射到拓扑上查看。

同样这两个模式，在全量规模上照样运行。随附的 preprint 示例使用了 302 种哺乳类、2,275 棵 free 拓扑与 2,275 棵 fix 拓扑基因树：

在共享的一致性子问题上，SplitAligner 的 Support(b) 与 PhyParts 高度一致（299 根内部枝上 Pearson r = 0.99997）。在论文报告的同机测试中，SplitAligner 完整流程用时 9 分 52 秒，PhyParts 的快速一致性设置用时 1 小时 50 分 53 秒——在这套数据上约快 11 倍。这是校准与运行时间比较，并不表示两种方法功能等价。

这里只让类群缺失发生，不引入基因树不一致。完整的 B1–B17 坐标轴始终冻结；每次删除后，由图决定的台账会更新为 observed、NA_fuse 或 NA_struct。在这个实验里，NA_topo 始终为零。