机器学习赋能无对齐k-mer频率距离计算与系统发育定位新方法kf2vec

本文推荐一种基于机器学习的新型系统发育定位工具kf2vec，它利用k-mer频率向量通过深度神经网络学习与参考系统发育树路径长度匹配的距离度量（Φ:k-mer→?m），无需序列比对（alignment-free）和标记基因，支持长读长、contig和基因组等多样数据类型的高精度放置（placement error 1.63 edges），在微生物和真核数据集上均优于传统方法（如CAFE的D2*和Skmer），为宏基因组分类和进化分析提供高效解决方案。

　　

1 引言

系统发育重建不仅用于理解进化历史，还能识别生物样本组成以解析生态环境。样本可能是未知个体（Bohmann等2020）或混合物（如微生物组研究），需通过已知物种序列数据库进行表征，包括分类学命名和系统发育分析（Washburne等2018）。尤其当库中缺少样本近缘物种时，系统发育定位（phylogenetic placement）——将新序列添加到现有参考树（Matsen等2010；Berger等2011）——成为解决方案，并可扩展至大数据集（Balaban等2022；Wedell等2023）。尽管定位方法众多（如Matsen等2010；Barbera等2019），但多针对短对齐序列（如读长或标记基因），且预处理复杂，涉及同源基因查找、对齐和组合结果（Mai和Mirarab 2022）。

无组装和无对齐定位成为有吸引力的替代方案，已有方法如基于k-mer的距离（Balaban等2020）或间隔词匹配（Lau等2019；Blanke和Morgenstern 2020）显示可行性，但精度仍不及最大似然法（Balaban等2022；Hasan等2022）。这些方法简化分析流程，并有望统一不同数据类型（如长读长、基因组、contig）到同一棵树。本文目标正是实现长输入序列的标记自由定位。

基因组距离与进化历史相关，距离-based系统发育推断有悠久历史（Warnow 2017），虽不如最大似然精确（H?hl和Ragan 2007；Bogusz和Whelan 2016），但继续应用，尤其距离-based定位更具扩展性和通用性（Balaban等2022）。计算系统发育意义距离的方法备受关注（如Tang等2019；Sarmashghi等2019；Ondov等2016；Shaw和Yu 2023；Sapci和Mirarab 2025）。

许多无对齐距离方法基于k-mer及其相关概念（Zielezinski等2017；Ren等2018），如长k-mer（k>20）存在/缺失或短k-mer（k<10）频率。方法差异在于将k-mer统计转化为基因组距离的公式（Zielezinski等2019），其理论基础和实证精度受关注（Vinga和Almeida 2003；Reinert等2009）。长k-mer方法包括Skmer（Sarmashghi等2019）、AFNN（Tang等2019）和Mash（Ondov等2016）；短k-mer频率有CAFE软件包实现28种无对齐相异度度量（Lu, Tang,等2017），其中背景调整k-mer计数的相异度最接近树距离，可用于下游系统发育分析（Lu等2021；Yuval Bussi和Ziv Reich 2021），但这些度量未被广泛采用，因其与给定树的相关性有限。

本文采用不同方式：不预定义公式，而是训练机器学习模型从参考数据集学习距离计算。提出的kf2vec（k-mer frequency to vectors）将输入（基因组、scaffold、基因组skim或长读长）表示为规范k-mer频率向量（维度为4^k/2），利用参考树学习将k-mer频率转化为嵌入，使嵌入间欧氏距离平方匹配参考树路径长度。新查询通过模型计算与参考序列的距离，用于距离-based定位（如APPLES-II）或分类分配。kf2vec并非首个从k-mer频率学习序列嵌入的机器学习方法，但与DEPP（Jiang, Balaban,等2022）有概念相似性。不同在于，kf2vec训练后需最少预处理，对参考树来源不可知，适用于多种测序数据类型。验证实验显示kf2vec在距离计算和系统发育定位上优于现有k-mer方法。

2 材料与方法

2.1 kf2vec：映射k-mer频率到系统发育距离

2.1.1 学习距离

设K=4^k/2（小k，如k=7），为所有规范k-mer数。任何足够长序列（组装、未组装短读集或单长读长）通过JellyFish（Mar?ais和Kingsford 2011）提取k-mer频率，归一化计数至和为1，因此每个输入表示为K维单形体向量（0≤v_i≤1, ∑₁^Kv_i=1）。计算问题：给定n个参考基因组，表示为n×K矩阵S，每行S_i,.为序列i的归一化k-mer频率向量。另有参考系统发育树T，叶标记为n个参考类群，任意两叶i和j的支系距离为T_ij（图1A）。固定嵌入维度m，寻求函数Φ:Δ^K-1→(?⁺)^m，使得∥Φ(S_i,.)-Φ(S_j,.)∥₂匹配√T_ij（1≤i,j≤n）。此目标类似DEPP（Jiang, Balaban,等2022），但输入为k-mer频率（无对齐），而非DEPP所用对齐序列，避免了多序列对齐（MSA）和同源标记基因检测，极大减少预处理时间。

函数Φ仅当可泛化时有用：若部分n序列在构建Φ时隐藏，∥Φ(S_i,.)-Φ(S_j,.)∥₂仍应匹配√T_ij。训练Φ后，用于未见查询序列计算其与参考树距离。对任何查询q，可计算n维距离向量d^(q)，包含d_i^(q)=∥Φ(q)-Φ(S_i,.)∥₂²值，给出查询到骨干序列的距离。

使用√T_ij而非T_ij，因de Vienne等（2012）和Layer和Rhodes（2017）显示对任何树T，存在点集P={Φ(t_i)}_i=1ⁿ在?^n-1欧氏空间中，点间距离对应√T_ij。因此对m=n-1，?₂范数合理；当m<n-1时，距离会有失真，但若失真小于分支长度则可容忍。

用机器学习求解，定义均方误差（MSE）损失：

Loss(Φ; S)=∑_i=1ⁿ∑_j=1ⁿ(∥Φ(S_i,.)-Φ(S_j,.)∥₂-√T_ij)² (1)

但长距离难估计且在距离-based系统发育中常降权（Fitch和Margoliash 1967；Beyer等1974），且对定位不重要（Balaban等2022）。因此仿Jiang, Balaban,等（2022），在MSE损失中降权长距离：

Loss(Φ; S)=∑_i=1ⁿ∑_j=1ⁿ1/(T_ij+p)(∥Φ(S_i,.)-Φ(S_j,.)∥₂-√T_ij)² (2)

其中1/T_ij为受Fitch和Margoliash（1967）启发的权重，使每项无标度，p为伪计数避免T_ij=0时权重无限。

构建Φ使用神经网络：两个前馈线性层，间隔ReLU激活函数。首层隐藏维度2m，次层m（默认m=1024）。末层输出为嵌入。用Adam优化器（Kingma和Ba 2014）优化损失函数，批量大小默认16（导致120对），运行固定轮数后选训练误差最低模型。学习率初始0.000013，逐步下调但不低于0.000003。用PyTorch（Paszke等2019）实现。默认kf2vec参数：m=1024，k=7（8192规范k-mer）。

2.1.2 用距离进行样本识别

2.1.2.1 系统发育定位

对每个查询q，估计距离向量d^(q)作为输入到距离-based定位方法APPLES-II（Balaban等2022），在树T上找最优放置（图1B）。APPLES-II用动态规划按最小二乘准则找q在T上最优位置，降权大距离。使用设置-f 0 -b 5，忽略长距离并用每个查询与骨干间最小5个距离。

2.1.2.2 分类分配

开发启发式算法基于嵌入距离分配分类标签。算法寻找与查询距离统计不可区分参考子集（称为邻居）。为检验不可区分性，假设到真实距离d的参考基因组估计距离乘以基因组长度L为泊松分布绘制，λ=d·L，即d?~Poisson(d·L)/L。对两参考估计距离d?₀≤d?_i，可计算p值：零假设最小为真值但较大为同分布噪声绘制，1-F_Poisson(d?_i·L; λ=d?₀·L)。基于此模型和额外准则决策：按距离升序迭代参考基因组，比较每个d?_i与最接近匹配d?₀和前一个d?_i-1，当满足两准则之一时停止并选{0,…,i-1}：要么比较d?_i与d?₀的p值低于0.1（即d?_i统计显著高于d?₀），要么观察到距离梯度远高于前迭代：d?_i-1-d?_i>3Δ_i，其中Δ_i=max_j=1^i-1(d?_j-1-d?_j)为至今观测最大增加（定义Δ₁=0.05×d?₁因无先验梯度）。

分类标签通过分析选中邻居中不同标签出现频率分配。从种水平开始并进展到更高级分类阶元（必要时），如果分类标签在至少3/4邻居中出现则分类到该标签。因此仅当较低阶元未达3/4共识时过程才上移到更高级分类。若无阶元满足阈值，查询分类到根类。

2.1.3 分治实现扩展到大树

构建单一模型面临两难：树表示为n×n矩阵，表大小二次增长，当n?10⁴时难存内存和优化神经网络过所有(ⁿ₂)条目。且随数据集系统发育多样性增加，单一模型可能难捕获所有复杂性。Jiang等（2023）提出用分治解类似问题：将树划分为较小子树，为每个子树训练模型，并用分类器为每个查询选最佳子集。采用类似方法：用TreeCluster（Balaban等2019）划分树为目标相似大小簇（用户指定最大，图1C）。然后为每个子树1…C分别训练嵌入函数Φ₁…Φ_C。还需训练分类器Ψ:Δ^K-1→Δ^C-1概率分配k-mer频率向量到子树。用标准交叉熵损失函数训练分类器。

查询时，先应用分类器Ψ选最高概率子树模型（图1C）。用该模型计算该子树所有参考序列与查询间距离。其他参考基因组距离未定义，APPLES-II允许。分类器模型Ψ架构与Φ共享首层，后接分类层（线性层加softmax函数产生概率）。

当n足够小，构建单一模型可行；但当数据集足够多样以致单一模型可能无法捕获全部复杂性时，用户可选分治，称为“包覆”模型（注意任何子树可用，不必是clade）。默认仅对n>4000参考数据集执行分治。

2.1.4 处理不完全数据（如contig、长读长）用分块

训练模型于全基因组计算k-mer频率对不完全数据（如长读长或contig）嵌入和放置非最优。解决方案涉及生成基因组许多子序列并训练模型于这些子序列。此设置让模型观察原始基因组变长区域并学习映射这些子序列到嵌入空间同一位置。为高效实现，用分块方法近似生成子序列。

首先将每个骨干基因组分成一组片段（块）以覆盖序列（图1D），默认块大小R=10,000 bp。块大多不重叠，但让块有短重叠以适应不被R整除基因组大小。用seqkit（Wei等2016）执行此步：先计算每个contig长度，再计算适当重叠使滑动窗均匀覆盖序列，最后提取每个结果块。排除总长<50 kbp基因组因高度不完全。对基因组i的每个块j，预计算其未归一化k-mer频率S′_i,j∈?^K，表示每个骨干基因组为一系列N有序k-mer频率向量，S′_i,1…S′_i,N。用此样本表示为输入训练分块版分类器和嵌入器。

训练中，每轮对每个基因组i，随机选连续子集k-mer频率向量S′_i,j…S′_i,j+l-1，首块j均匀随机选。采样向量数l从指数分布绘制，均值设N/5（20%总块数），截断在[1,N-j)内。选择截断指数分布确保高方差（如若N=100和j=50，μ≈16.4，σ²≈400）。采样子集然后作为连续块l向量始于j在块集内。然后选中块k-mer频率相加并重新归一化创建单个输入项S_i.,=1/(l·(R-k+1))∑_a=0^l-1S′_i,j+a。此过程每构建批时重复，使每轮基于新子序列原始基因组训练。此分块策略用于训练分类器和嵌入器。对嵌入器（非分类器），为鼓励模型映射不同子序列到同一嵌入，每批对每个基因组包含两轮子采样（即S′_i,.和S″_i,.用不同j和l构建），设相同基因组两子样本间距离为0。结果分块模型批大小比未分块大两倍，导致约四倍多对（496对120，默认）。当包覆和分块同用，所有相同基因组块训练嵌入器时分配到相同真clade。

推断时无分块。每个输入基因组、contig或长读长整体用作方法输入，计算一个k-mer频率向量。因查询时输入未分块，整个输入分类到一个clade并嵌入到一个位置。因此包覆和分块可同时使用。

2.2 评估实验

在九个不同属性数据集（表1）上系列实验评估提出方法，聚焦几个问题。每个案例有参考树。多案例执行留出实验，其中树叶从训练隐藏并用作查询。这些情况，视查询在参考上放置为地面真值。作为误差度量，比较方法所得放置与地面真值，报告其间节点数，称为放置误差。有些案例查询放置在树T₁与参考树T₂不同。这些案例，比较更新树与T₂的Robinson和Foulds（RF）距离（Robinson和Foulds 1981）（两者诱导有骨干和一查询），减去T₁和T2间RF距离，报告结果称为delta误差；delta误差显示放置产生额外树距离，当T₁=T₂时匹配放置误差。接下来描述每个实验，细节和确切命令留附录S3。

2.2.1 D1/D2：全长微生物基因组

为在微生物数据基准测试工具，用Web of Life数据集，含10,575类群系统发育树（Zhu等2019）（WoL19）和新版199,330基因组（Balaban等2023）（WoL23），均从约380基因树构建。以两种方式选未见训练数据查询（D1和D2在表1）。初始分析聚焦检查kf2vec算法变体，从WoL19随机选500查询集移除并用余下为参考（D1）。在比较kf2vec与其他方法分析中，用全WoL19树为骨干树并选500出现在WoL19但缺WoL23查询为查询（D2），报告delta误差。这些查询选有范围距离到WoL19最近参考类群，允许检查新颖性影响。定义新颖性为查询到参考数据集中任何基因组最小基因组核苷酸距离，用Mash（Ondov等2016）近似。

2.2.1.1 kf2vec参数影响

默认参数对D1数据集表现好（图2B）。探索替代k值，观察k=7在原核生物测试得最低放置误差（图2B）。更短和更长k-mer导致更高放置误差。虽k=8近7且可能对其他基因组（如更大真核基因组）提供更好精度，选k=7为默认因更低运行时间和更高精度。加一或两个额外线性层带非线性激活不改进精度但显著增加训练时间40%和130%（图S1B）。

2.2.1.2 训练

接下来检查训练过程参数。训练和测试损失首100轮急剧下降（图S3A）。训练损失继续降8000轮。但测试损失超约1000轮后无改进。尽管如此，稳定训练损失下，放置误差继续改进至少4000轮，最实质增益在首1000轮（图S3B）。给定此趋势，建议训练4000-8000轮足够。增加学习率和/或最小学习率10×有小但负影响放置误差，而100×增加大幅恶化结果（图S1B），显示选适当学习率重要。还测试增加学习率衰减和批大小2×会否影响模型收敛，但这些变化无测量效应放置误差。

2.2.1.3 分治（包覆）

包覆有效减少误差和运行时间。总体，在子树训练嵌入器模型优于未包覆模型（平均误差1.81 vs. 1.63），无论查询新颖性（图2A）。用于分配到子树分类器对98.2%查询准确。错误分类案例 among最 novel查询，距离最近基因组范围0.21到1.6（平均0.98）。若放置每个查询在其真子树以聚焦嵌入器，见包覆模型在15子树中13个得更好放置（图S4A）。放置误差跨子树变化，范围1.1到2.5边，显示有些树部分更具挑战性。显著，误差差异 between clades不总由clade大小或子树直径解释（表S4）；较高平均误差clades有正常中值误差，显示其较高平均误差 mostly due异常值。除精度轻微改进，包覆还通过允许clades独立训练帮助可扩展性，将看到。

2.2.1.4 过滤参考

基于查询质量度量影响，可能怀疑参考基因组质量过滤也可能改进精度。测试此想法并观察全基因组放置精度保持不变即使模型