《Cladistics》：Preprocessing missing data and partitions for phylogenetic analysis under dynamic homology, with application to museomics and phylogenomics
编辑推荐：
在动态同源性分析中，数据通常必须经过预处理，以区分测序伪迹（缺失数据）与插入–缺失事件，但使用POY/PhyG的实证研究此前多采用人工预处理方式。本研究比较了动态同源性分析中处理缺失数据的不同流程，提供了实现这些流程自动化的Python脚本，并使用60个模拟数
在动态同源性分析中，数据通常必须经过预处理，以区分测序伪迹（缺失数据）与插入–缺失事件，但使用POY/PhyG的实证研究此前多采用人工预处理方式。本研究比较了动态同源性分析中处理缺失数据的不同流程，提供了实现这些流程自动化的Python脚本，并使用60个模拟数据集和4个实证数据集（包括Sanger测序、系统发育基因组学及历史DNA序列）评估了其性能。研究人员阐明了将缺失数据编码为间隙（gap）或未知核苷酸/氨基酸在理论和方法学上均是不恰当的。相比之下，通过分区方式隔离缺失数据，或将其编码为任意观测状态或无状态（“?”，即最大多态状态；核苷酸为N+间隙，氨基酸为X+间隙）是可行的，但存在权衡。分区可以降低时间复杂度，但将潜在同源核苷酸隔离在不同分区中，从而削弱了同源性的检验；反之，问号能够实现不受约束的同源性检验，但会增加时间复杂度。为平衡上述权衡，研究人员提出了选择最优分区的启发式策略。所述方法在大规模模拟数据集中降低了时间复杂度，并改善了所有实证数据集的树代价，凸显了在动态同源性分析中正确识别缺失数据的重要性。
**动态同源性分析中的缺失数据预处理与分区策略：方法开发与实证验证**
**一、研究背景与问题**
系统发育分析中，序列数据的同源字符与字符状态识别存在两种基本范式：静态同源性（static homology）与动态同源性（dynamic homology）。静态同源性先通过多序列比对（multiple sequence alignment, MSA）基于与树最优性解耦的相似性函数确定同源字符，再在树上优化该比对并搜索最优树。动态同源性则在搜索最优树的过程中直接将未比对序列优化到每棵候选树上（即树比对，tree alignment），同源字符与状态在事后通过最优树各节点的转化列表或隐含比对（implied alignment）确定。在组合优化中，给定未比对序列集合、转化代价与最优性准则下寻找最优树的问题被称为广义树比对问题（Generalized Tree Alignment Problem, GTAP），其与TAP均为NP-hard问题，精确解对大多数实证数据集不可行，因此已有直接优化（direct optimization, DO）、迭代通路优化与基于搜索的优化等启发式算法被提出。
实证研究中，数据缺失可能源自演化上的插入–缺失事件（indels），也可能源于技术限制。由于树比对将所有序列长度变异视为真实并经indels解释，故必须在树比对前准确识别并编码伪迹缺失数据，以免终端类群因测序、组装和拼接伪迹而错误聚类。例如，引物或探针失效、拷贝数变异等可导致整个位点缺失；侧翼低质量碱基剔除、不同引物组无重叠扩增、高通量测序覆盖不足等可导致位点部分缺失。尽管输入序列在树比对下未比对，多数实证分析仍借助初步相似性比对来标识技术性缺失数据。当各位点单独输入时，缺失位点对相应终端被处理为零长度分区。然而，位点内部长度变异需通过三种途径预处理并各有取舍：其一，裁剪完整序列以匹配不完整序列长度，但会排除潜在重要证据；其二，延伸不完整序列以匹配完整序列长度，早期POY版本通过插入未知身份核苷酸（IUPAC N）实现，但该方式假设缺失序列与完整序列等长，未必正确，POY 5（Wheeler et al., 2015）与PhyG（Wheeler et al., 2024）遂允许以问号（?）扩展，表示任意观测状态或无状态（即最大多态状态；核苷酸为N+或间隙，氨基酸为X+或间隙）；其三，以井号（#）隔离缺失数据区块，将序列打断为更小片段，可加速分析并降低内存需求，但会限制相邻片段间同源性检验，可能导致次优解。已有研究多为人工预处理，未明确评估相关选择对同源性、树代价或拓扑结构的影响，致使动态同源性中缺失数据处理的最佳途径尚未明确；人工预处理亦增加主观性、降低可重复性且易出错。在包含成百上千基因比对的系统发育基因组学以及具有高比例变异分布缺失数据的历史DNA（hDNA）和古DNA（aDNA）研究中，实现自动化预处理尤为必要。
**二、研究内容**
针对上述问题，研究人员（1）比较了树比对前处理缺失数据的多种流程，（2）提供了自动化输入序列预处理的Python脚本（prepDyn），（3）利用模拟和实证数据及DO算法对该脚本进行了测试。理论上，将缺失数据编码为IUPAC N不恰当地假设了未知核苷酸实际存在且未考虑长度变异可能；以问号编码缺失数据则不作该假设，可减少ad hoc假设，对hDNA序列尤为重要。问号与分区均可恰当编码缺失数据，但分区将同源检验限制在同一分区内，比问号施加了更多假设。分区虽增加约束，却可降低时间复杂度，使更彻底的树搜索成为可能。为平衡二者，研究人员提出四种启发式分区方法：最大分区（-pm max，在所有缺失数据区块两侧添加井号列）、平衡分区（-pm balanced，基于逐轮分区的合并阈值）、保守分区（-pm conservative，在高度保守区域插入井号）和等长分区（-pm equal，按指定数量将序列划分为等长片段），并以逐轮分区策略确定最优分区数。
**三、关键技术方法**
研究人员开发了整合四步流程的prepDyn.py脚本：数据收集、裁剪、缺失数据识别和逐轮分区。数据收集支持从GenBank获取序列或解析本地多序列比对文件，并可组合多来源数据，初步比对使用MAFFT或Clustal W生成。裁剪包括去除侧翼不变列和“孤儿核苷酸”（被大量间隙隔离的短侧翼核苷酸块，阈值x通常为25–45）。缺失数据识别将末端间隙自动替换为问号，非重叠扩增片段间的间隙转换为问号，并可将hDNA中的IUPAC N替换为问号。分区步骤实施上述四种启发式方法。模拟研究中，研究人员使用ngesh模拟了含10、20、40、80个叶子节点的系统发育树，以iqtree3 --alisim按JC模型模拟了100、1000和10000个核苷酸的DNA序列，以等长分区控制分区数量，在PhyG v.1.3中以build和swap命令执行DO分析并进行回归分析；实证研究使用了4个数据集：新热带湍蛙科（Hylodidae，75叶子，约4291 bp，4个Sanger片段，Sá et al., 2022）、羊膜动物（Amniota，16叶子，约187026 bp，248个蛋白编码基因，Chiari et al., 2012）、新热带树蛙（Dendropsophus，115叶子，约6304 bp，8个片段含Sanger和hDNA序列，Nakamura et al., 2025）以及鲈形鱼类（Percamorpha，14叶子，约99948 bp，100个锚定杂交富集片段，Dornburg et al., 2017）。每数据集采用六种策略（无预处理、预处理不分区、预处理加最大/平衡/保守/等长分区），共执行252次等权简约分析，组合使用POY v.5.1.1、PhyG v.1.3和TNT v.1.6，并计算归一化Robinson-Foulds距离和分支长度差异比例。
**四、主要研究结果**
*模拟分析*：prepDyn预处理运行时间为0.01–2.07 CPU秒，回归模型（校正R2=0.99）显示运行时间随叶子数与总序列长度近似线性增长（估计值分别为0.88和0.95，均P<0.001），分区数非显著预测因子（P=0.37）。DO下游分析运行时间为188–42494920 CPU秒，时间复杂性符合O(n^3–4m^x)预期。在最大模拟数据集（80叶子，10000核苷酸）中，分区数轻微降低运行时间（估计值=–0.16，P=0.002）；在最小数据集中分区数轻微增加运行时间（估计值=0.25，P=0.005），即分区在大数据集中可缩短运行时间，但在小数据集中不然。树代价与分区数呈显著正相关（估计值=1.78×10^–3，P=0.02）。
*实证分析*：所有经预处理的策略（策略2–6）树代价均低于无预处理策略（策略1）。最优分区策略因数据集而异：Dendropsophus数据集中以等长分区最优（13 461步），Amniota数据集中亦以等长分区最优（194 372步，第4轮），Hylodidae数据集中以策略5保守分区最优（4 746步），Percamorpha数据集中以策略2不分区分区最优（70 443步）。逐轮分区在降低树代价方面发挥重要作用。归一化RF距离显示，系统发育基因组学数据集（Amniota和Percamorpha）中策略2–6间多数树对拓扑完全相同（RF=0），但Sanger数据集（Hylodidae和Dendropsophus）中拓扑存在差异；分支长度差异则相反，系统发育基因组学数据集（77–100%）高于Sanger数据集（18–78%）。
**五、讨论与结论**
讨论部分指出，prepDyn自动化的预处理流程对大规模系统发育基因组学数据及hDNA序列尤为重要，能够避免人工处理的繁琐与易错性。已有工具如Gblocks和trimAl仅侧重裁剪，缺乏动态同源性所需的分区功能以及对非标准符号（问号和井号）的处理能力。研究人员澄清了POY/PhyG中“零长度分区被解释为缺失数据”并非缺陷，而正是程序的设计意图；用户需告知程序哪些变异源于技术伪迹。理论上，基因内分区应最小化以增强解释力，但分区可加速实证分析并扩展树比对空间的探索，在时限搜索中尤为有价值。保守、等长和最大分区在不同数据集中各有优势：保守分区对同源检验约束最小，等长分区通过改变井号位置多样化搜索空间，最大分区则最大程度降低计算负担。值得注意的是，尽管两个系统发育基因组学数据集在各分区策略下拓扑相同，但大多数或全部分支长度存在差异，揭示分区策略对同源性的影响在拓扑不变时仍可影响依赖分支长度的下游分析（如多样化速率估计、系统发育回归）。
研究结论可归纳为三点：第一，prepDyn以问号（而非IUPAC N和间隙）为主要编码方式实现DNA序列预处理自动化；若执行分区，缺失数据亦可编码为零长度分区；若跳过分区，prepDyn可推广至将间隙视为第五字符状态的其他系统发育程序。第二，设计系统发育分析策略时，除叶子数和核苷酸数外，分区数亦应纳入时间复杂度的预测变量；实证数据集中需多次限时搜索，逐轮分区可改善树代价，prepDyn提供的四种分区方法（平衡、保守、等长、最大）可通过树代价进行比较。第三，prepDyn将先前人工执行的步骤自动化，降低了主观性和人为错误，提高了可重复性。
