PhageCGRNet：将基因组的混沌游戏表示与卷积神经网络相结合，以实现精确的噬菌体宿主分类预测

《PLOS Computational Biology》：PhageCGRNet: Integrating Chaos Game Representation of Genomes with Convolutional Neural Network for accurate phage host classification prediction

　　摘要：噬菌体（或细菌病毒）在微生物群落中扮演着关键角色，准确预测噬菌体的宿主对于理解这些病毒的动态及其对细菌种群的影响至关重要。在预测噬菌体宿主分类的过程中，特征提取是一个关键步骤，它直接影响预测的准确性。在用于特征提取的技术中，k-mer是最常用的方法。尽管已经提出了许多基于k-mer的方法，但这些方法通常仅使用k-mer的频率信息作为特征。然而，当频率相同时，这些k-mer的频率信息就变得不那么有用了。为了解决这一限制，我们提出了一种名为PhageCGRNet的新方法，该方法不仅利用了k-mer的频率信息，还结合了k-mer的位置信息。在我们的方法中，我们将每个基因组序列表示为一个三维矩阵，其中包含k-mer的频率特征和位置特征，然后使用卷积神经网络模型来预测宿主类别。具体来说，我们将直接从序列中提取的k-mer的频率信息与使用混沌游戏表示方法获得的k-mer的位置信息结合起来构建特征矩阵，该矩阵作为卷积神经网络的输入。我们在两个基准数据集上进行了实验，并将PhageCGRNet与现有的先进噬菌体宿主分类方法进行了比较。实验结果表明，与其他最先进的方法相比，PhageCGRNet在这两个数据集的物种和属两个分类层次上都取得了更高的准确性。

作者总结：随着耐药性细菌病原体对主要临床使用抗生素的普遍性增加以及抗生素研究和开发进展缓慢，噬菌体能够特异性地靶向并杀死细菌的能力使其在抗菌治疗中变得越来越重要。为了实现有效的噬菌体疗法，准确预测噬菌体与其宿主之间的相互作用是必不可少的。然而，大多数现有的预测方法仅依赖于k-mer的频率信息，常常忽视了序列中位置信息的重要性。为了解决这一限制，我们在本研究中提出了一种新的噬菌体宿主预测方法。该模型在两个基准数据集上的宿主分类性能优于当前最先进的方法，为推进噬菌体疗法的实际应用提供了强大的计算工具支持。

1. 引言：噬菌体被认为是地球上最丰富且可能最多样化的微生物[1]。噬菌体通常具有特定的宿主，它们会结合到宿主细胞表面的特异性受体上，并将遗传物质注入宿主细胞[2]。由于噬菌体能够精确高效地攻击宿主细菌而不伤害人类或其他生物，因此可以用来选择性地杀死病原菌，从而克服抗生素耐药性带来的挑战[3]。同时，通过改变噬菌体的DNA序列，可以设计出用于基因传递、基因治疗或肿瘤治疗的定制噬菌体[4]。因此，研究噬菌体非常重要，但传统的实验室培养噬菌体需要大量的时间和资金[5]。此外，不到1%的宿主已在实验室中被培养出来[6]。下一代测序技术为发现新的噬菌体提供了新的方法[7]。噬菌体测序数据的急剧增加使得宿主鉴定变得困难，并限制了噬菌体宿主分类的研究。因此，有必要开发更准确的噬菌体宿主分类工具。近年来，提出了多种计算方法来预测噬菌体宿主分类，这些方法大致可以分为基于比对的方法和不基于比对的方法。基于比对的方法可以通过两种主要策略实现。第一种策略涉及识别噬菌体和宿主之间的共享片段[8]。这些方法包括匹配CRISPR间隔序列、识别宿主基因组中的相关噬菌体、噬菌体和宿主之间水平转移的基因。然而，使用间隔序列进行噬菌体宿主预测存在局限性，而且许多宿主没有间隔序列[9]。第二种策略依赖于噬菌体之间的相似性[10]。这些方法不是直接寻找噬菌体-宿主匹配，而是通过将查询噬菌体与已知宿主的噬菌体数据库进行比较来推断宿主，基于序列相似性。例如，RaFAH[11]通过使用隐马尔可夫模型搜索噬菌体蛋白与蛋白质家族数据库，并将随机森林分类器应用于得到的比对分数来预测宿主。类似地，vHULK[12]也依赖于蛋白质相似性进行宿主预测。它将噬菌体蛋白与原核病毒数据库进行比较。比对显著性分数被转换并用作输入特征。最后，使用多层感知器神经网络来预测宿主。不基于比对的方法包括丰度分析方法和核苷酸组成方法[13]。其中，核苷酸组成方法基于这样的观察：噬菌体及其宿主在密码子使用或k-mer（k-mer表示长度为k的DNA片段）上经常表现出相似的模式。此外，感染同一宿主的噬菌体由于宿主特异性选择或共同的进化历史而倾向于共享相似的k-mer特征。因此，k-mer频率的相似性在噬菌体宿主预测中被广泛使用。这些k-mer可以作为宿主分类算法中的特征提取技术。提取的特征随后被输入到机器学习算法中以进行噬菌体宿主分类[14–21]。一种仅依赖于噬菌体序列的突出方法是DeepHost[22]。它将序列编码为3D矩阵，并使用卷积神经网络模型来分类噬菌体的宿主。相比之下，其他几种方法同时利用了噬菌体和宿主的数据。例如，CHERRY和PHPGCA在噬菌体和宿主之间建立了网络。具体来说，CHERRY[23]构建了噬菌体和宿主之间的知识图谱，然后使用图卷积神经网络模型来预测宿主的物种。PHPGCA[24]也在噬菌体和宿主之间构建了异构图谱，然后使用图对比学习算法增强节点表示，将宿主预测任务视为推荐任务。另一种方法CL4PHI[25]使用频率混沌游戏表示（FCGR）[26]来提取序列特征，然后使用对比学习算法来识别感染同一宿主的噬菌体以预测噬菌体宿主关系。不同地，PB-LKS[27]利用噬菌体和宿主之间最相似的基因组片段的局部k-mer相似性作为特征。然后使用随机森林模型来预测噬菌体与其宿主之间的相互作用。虽然一些现有方法，如DeepHost[22]和CL4PHI[25]，已经展示了将k-mer频率图与卷积神经网络（CNN）模型[28]结合用于宿主预测的实用性，但它们主要依赖于频率信息。相比之下，我们方法PhageCGRNet的关键创新在于其新颖的特征表示。我们将k-mer频率与来自混沌游戏表示（CGR）的位置信息整合到一个统一的三维特征矩阵中。这种整合使我们的模型能够同时学习哪些k-mer存在以及它们在基因组序列中的位置，从而捕获更全面的序列特征集。在这项工作中，我们提出了一种新的噬菌体宿主分类方法，称为PhageCGRNet。除了直接从基因组序列中提取特征外，我们还使用CGR[29]来重新提取特征。具体来说，我们通过考虑序列中k-mer的频率和k-mer在CGR图中的相应平均位置来获取序列信息。然后将这些整合的信息作为特征矩阵输入到CNN模型中，该模型学习噬菌体的特征。这些学到的噬菌体特征随后被用来预测它们对应的宿主。特征矩阵中k-mer的频率和位置信息的结合提供了更全面的表示，然后与深度学习算法CNN结合，以更准确地进行预测。值得注意的是，我们的方法使用k-mer频率和位置信息来表征噬菌体基因组。宿主预测是通过比较噬菌体之间的序列差异来识别与特定宿主相关的基因组特征，而不是通过噬菌体与其宿主之间的直接序列比对来实现的。

2. 结果：2.1. 数据集：在这项工作中，我们使用了DeepHost[22]和CHERRY[23]之前使用过的两个数据集，分别称为数据集1和数据集2。数据集1最初包含8,756个噬菌体。这些噬菌体根据DeepHost[22]中描述的数据选择方法进行了过滤。过滤后，我们保留了8,595个属于49个属级和72个属级宿主分类的噬菌体序列，以及数据集1中属于118个种级宿主分类的7,483个噬菌体序列。数据集2包含1,875个噬菌体，涵盖了206个种级、108个属级和67个属级宿主分类。对于数据集1和数据集2，我们以8:1:1的比例随机划分它们，以获得训练集、验证集和测试集。2.2. 性能评估指标：为了将我们的模型与最先进的方法进行比较，我们使用准确性、精确度、召回率、F1分数和top-k准确性来衡量模型的性能。该准确性指标的公式在方程（1）中给出：(1)宿主分类是一个多类预测问题。因此，我们使用宏观精确度、宏观召回率和宏观F1分数来进行性能评估。计算公式如下：(2)其中，是实际属于类别i并被正确预测为类别i的样本数量。是实际上不属于类别i但被错误预测为类别i的样本数量。(3)其中，是实际属于类别i但被错误预测为另一个类别的样本数量。(4)2.3. 性能比较：2.3.1. 基线方法：我们在两个基准数据集上进行了实验：数据集1（Cherry）和数据集2（DeepHost）。为了全面评估我们方法的性能，我们将其与几种最先进的方法进行了比较。这些方法包括RaFAH[11]和vHULK[12]、DeepHost[22]、CHERRY[23]、PHPGCA[24]、CL4PHI[25]和PB-LKS[27]。对于RaFAH[11]、vHULK[12]、PHPGCA[23]和PB-LKS[26]，我们使用原始论文中提供的参数设置重新进行了数据集1上的实验。对于其他三种方法（DeepHost[22]、CHERRY[23]和CL4PHI[25]），我们使用了CL4PHI[25]在数据集1上报告的结果。对于数据集2，由于数据集划分的变化，我们使用它们各自论文中指定的参数设置重新实现了所有七种方法。2.3.2. 实验设置：在我们的实验中，特征提取阶段使用的k-mer长度设置为7。这导致了一个大小为3×128×128的特征矩阵。我们的CNN模型由两个卷积层组成。每个卷积层的核大小为5×5，步长为1，填充为2。每个卷积层之后是一个核大小为2×2、步长为2的池化层。此外，全连接层包含512个神经元。我们在表1中提供了详细的参数设置。同时，我们将我们提出的PhageCGRNet与其他先进方法进行了比较，包括RaFAH[11]、vHULK[12]、DeepHost[22]、CHERRY[23]、PHPGCA[24]、CL4PHI[25]、PB-LKS[27]。我们使用四个指标（准确性、精确度、召回率和F1分数）在两个数据集的三个分类层次（物种、属和科）上全面评估了我们方法的性能。详细结果见表2。粗体值表示每个指标的最佳性能。由于RaFAH[11]在属级别进行预测，因此在物种级别用“-”表示。如表2所示，在数据集1上，我们的方法在所有指标上均优于所有其他最先进的方法。在物种级别，尽管精确度略低于CL4PHI，但我们的方法在其他所有指标上取得了更好的结果。在数据集2（表3）上，我们的方法在物种级别的准确率为79.68%，比第二好的方法PHPGCA高出6.42%；在属级别达到了89.30%的准确率，比PHPGCA高出8.55%；在科级别达到了91.44%的准确率，提高了4.76%。此外，我们的方法在其余评估指标上也始终优于所有其他方法。

下载：PPT PowerPoint幻灯片PNG更大图像TIFF原始图像表1. CNN架构的详细配置。https://doi.org/10.1371/journal.pcbi.1014130.t001下载：PPT PowerPoint幻灯片PNG更大图像TIFF原始图像表2。对数据集1上的宿主预测性能的全面评估以及与其他方法的比较。 https://doi.org/10.1371/journal.pcbi.1014130.t002 下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像表3. 对数据集2上的宿主预测性能的全面评估以及与其他方法的比较。 https://doi.org/10.1371/journal.pcbi.1014130.t003除了上述评估指标外，我们还使用了Top-k准确率（k=1,5,10,15,20）来评估预测性能。图1中的条形图显示我们的方法在所有k值上都取得了高性能。此外，图2中呈现的混淆矩阵详细可视化了模型在不同宿主类别间的分类行为。对角线元素代表每个类别的正确预测数量，而非对角线元素揭示了类别间的误分类分布。对角线元素颜色较深，表明大多数预测落在正确的类别上。非对角线元素稀疏且颜色较浅，证实了误分类率较低。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像图1. 数据集1和数据集2在物种、属和科级别的Top k准确率。 https://doi.org/10.1371/journal.pcbi.1014130.g001 下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像图2. PhageCGRNet的分类性能混淆矩阵。 https://doi.org/10.1371/journal.pcbi.1014130.g002 2.3.3. 统计测试。为了进一步证明我们的方法在数据集1上比CL4PHI方法更准确，我们评估了PhageCGRNet和CL4PHI的评估指标值之间的统计显著性。首先，我们分别重复实验30次。我们使用Kolmogorov-Smirnov检验在0.05的P值下测试正态性[30]。结果表明数据不符合正态分布。最后，我们采用了Wilcoxon符号秩检验[31]。PhageCGRNet的准确率显著高于CL4PHI（P值<0.05），且效应量较大（Cliff’s Delta [32] = 0.4775）。2.4. 在训练集和测试集相似性变化下的模型性能我们使用Dashing工具进行了分层相似性分析，以计算测试集和训练集之间的相似性。我们根据相似性（≤20%、≤40%、≤60%、≤80%和≤100%）对测试集进行了划分。然后在数据集2的物种分类级别内评估了所有竞争模型的预测准确率，结果如图3所示。分析显示，随着测试样本和训练集之间相似性的增加，所有模型的性能都有所提高，但我们的方法在每个相似性区间内都实现了最高的准确率。这证实了我们方法的优势是稳健的，并不严重依赖于训练集和测试集之间的特定相似性分布。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像图3. 训练集和测试集之间的相似性影响宿主预测性能。 https://doi.org/10.1371/journal.pcbi.1014130.g003 2.5. 消融实验为了确定影响结果准确性的因素，我们进行了消融实验。如表4所示。实验设置如下：首先，我们使用k-mer频率特征和CGR位置特征作为独立输入，通过多层感知器（MLP）模型[33]进行分类预测。随后，将这两个特征整合到一个复合的3D特征矩阵中，用于基于MLP的分类。同时，我们也使用卷积神经网络（CNN）重复了相同的工作流程，分别使用每个单独的特征及其组合形式（频率和位置特征）执行分类任务。结果表明：对于CNN和MLP模型，CGR衍生的位置特征在准确率上都有所提高（详见表4），与k-mer频率特征相比。此外，频率和位置信息的协同使用在进一步提高预测准确性方面发挥了不可或缺的作用。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像表4. 数据集1和数据集2中的消融实验。 https://doi.org/10.1371/journal.pcbi.1014130.t004为了探讨CNN是否有助于提高我们方法的准确性，我们进行了模型性能的比较分析，比较了CNN集成前后的性能（以MLP模型为基准）。实验框架包括两个维度：（i）单特征比较：使用k-mer频率特征矩阵和CGR位置特征矩阵独立评估CNN和MLP在单特征条件下的分类性能差异；（ii）复合特征验证：将两种特征类型融合到一个3D矩阵中，以验证混合输入下的模型行为。实验结果表明：在单特征场景下，CNN在两个数据集的物种级别/属级别/科级别分类任务中的表现优于MLP基线（具体数值见表4）。当使用组合特征时，CNN模型在捕捉特征交互方面表现出更强的能力，从而提高了两个数据集的准确率。这些结果明确表明，CNN集成有效地增强了特征表示学习，从而提高了分类准确性。2.6. 参数分析由于我们的方法基于k-mer，k的值是一个关键参数。为了确定其最优值，我们进行了实验研究。我们使用了一系列k-mer长度来测试我们的方法。对于每个k-mer长度，我们评估了相应的准确率。这些实验的结果显示在图4(A)中。从这个图中，我们可以清楚地观察到k值与预测准确率之间的关系。最初，随着k值的增加，噬菌体-宿主分类的准确率也随之提高。当k等于7时，准确率达到峰值。然而，当k继续增加超过7时，准确率开始下降。基于这些发现，我们得出结论，我们方法的最优k值为7。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像图4. 不同参数设置下的宿主预测准确率。(A) PhageCGRNet在数据集1和2上的宿主预测准确率。(B) 使用不同深度的卷积层在两个数据集上的宿主预测性能。在两个面板中，最优值用黑色标出。 https://doi.org/10.1371/journal.pcbi.1014130.g004此外，我们还研究了卷积层数量对模型性能的影响。为此，我们进行了不同数量的卷积层实验，以确定最优数量。如图4(B)所示，实验结果表明使用两个卷积层可以获得最佳性能。这可能是因为两个卷积层在模型复杂性和计算效率之间取得了良好的平衡，能够在不导致过拟合或过度消耗计算资源的情况下捕获足够的特征信息。此外，我们测试了四种不同的核苷酸到CGR顶点映射方案，以确定最优分配策略。由于核苷酸可以以多种方式分配给CGR顶点，我们特别检查了将（1, 1）点分配给不同核苷酸的影响。如表5所示，实验结果表明不同的分配方案对模型准确率的影响很小。因此，我们采用了以下映射方案：A（0, 0）、C（0, 1）、G（1, 0）和T（1, 1）。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像表5. 使用不同核苷酸分配方案的宿主预测准确率比较。 https://doi.org/10.1371/journal.pcbi.1014130.t005尽管不同分配方法之间的模型准确率存在轻微差异，但这些差异非常有限。为了合理解释这一现象，我们进行了一系列实验：每种分配方法独立重复了20次，每次实验使用相同的随机种子来控制变量。随后，我们使用配对t检验进行了统计显著性分析。如表6所示，我们发现所有分配方法的P值都显著大于0.05的显著性阈值，表明不同分配方法之间的模型性能差异没有达到统计显著性。在这项工作中，我们在后续实验中采用了ACTG分配方案作为默认配置。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像表6. 使用不同核苷酸分配方案的配对t检验的P值。 https://doi.org/10.1371/journal.pcbi.1014130.t006 2.7. CNN与其他机器学习模型的比较我们将PhageCGRNet的CNN模型与各种其他机器学习分类器进行了比较，包括传统机器学习领域的支持向量机（SVM）[34]和随机森林（RF）[35]，以及深度学习领域的循环神经网络（RNN）[36]。我们使用我们的特征提取方法为所有分类器获取输入矩阵。图5显示了每个模型的预测准确率比较。为了进一步定量评估CNN模型与其他分类器之间的性能差异，我们进行了统计显著性测试和效应量分析。基于30次重复实验的结果，我们使用Wilcoxon符号秩检验来评估统计显著性，并计算了Cliff’s Delta效应量（δ）来量化差异的大小。测试结果总结在表7中。从该表可以看出，在物种级别、属级别和科级别，P值均低于0.05，且多个Cliff’s Delta效应量达到或接近1。这表明CNN模型的准确率显著高于其他分类器，效应量也较大。从图5和表7中可以清楚地看到，CNN模型的准确率高于其他分类器。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像表7. CNN模型与其他分类器的统计比较。 https://doi.org/10.1371/journal.pcbi.1014130.t007 下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像图5. PhageCGRNet与其他机器学习分类器在数据集1和数据集2上的科级别、属级别和物种级别的噬菌体宿主分类性能比较。 https://doi.org/10.1371/journal.pcbi.1014130.g005 2.8. 模型可解释性和特征分析为了评估整合CGR位置信息后噬菌体宿主聚类的改进，我们对两种表示方法进行了定量聚类评估：一种仅依赖于k-mer频率特征，另一种结合了k-mer频率和CGR位置特征。具体来说，我们计算了它们的轮廓系数[37]和Calinski-Harabasz指数（CH）[38]。结果显示，融合的CGR和k-mer特征集的轮廓系数为0.2889，CH指数为126.55，这显著高于仅使用k-mer特征时获得的0.1216和59.44。这表明感染同一宿主的噬菌体在高维特征空间中形成了更紧密、更分离的簇，定量验证了整合CGR位置信息后的聚类效果增强。通过结合CGR的全局分布和k-mer的局部频率，我们有效地克服了单一特征的局限性。值得注意的是，它解决了k-mer本身无法表示重复元素的空间分布的问题，验证了融合这两种类型信息的优越性。Jeffrey [28]和Deschavanne等人[39]指出，CGR可以揭示基因组特有的全局序列模式，包括双勺结构、对角线分布、从上到下的强度变化或反向变化、空区域以及不同形状的富词区域。例如，在图6中，Gordonia、Pseudomonas和Mycolicibacterium的AT对角线区域的空白区域是由于缺乏连续的AA/TT碱基对（或以它们结尾的k-mer）和固定的组合AT/TA二核苷酸（或以它们结尾的k-mer）造成的。Staphylococcus和Flavobacterium中AT对角线的高密度着色反映了序列中周期性重复的AT/TA二核苷酸的存在。Flavobacterium的CG子象限中点密度的降低表明缺乏以CG结尾的子序列。这些模式揭示了序列的全局碱基组成偏好和分布特征。下载： PPTPowerPoint幻灯片PNG较大图像TIFF原始图像图6. 从数据集2中十个最丰富的噬菌体属生成的FCGR图。红色、青色和蓝色框分别突出显示了贡献最大的区域、第二大的区域和第三大的区域。 https://doi.org/10.1371/journal.pcbi.1014130.g006为了进一步解析PhageCGRNet模型识别的序列模式，我们使用了基于博弈论的SHapley Additive exPlanations（SHAP）值[40]来分解特征对模型预测结果的贡献。SHAP值不仅量化了每个输入特征对单个预测的边际影响，还确保了全局可解释性的一致性。具体来说，我们从数据集2中十个最丰富的噬菌体属中选择了基因组序列，计算了所有输入层特征的SHAP值，按贡献分数对这些特征进行了排名，并提取了每个属中贡献最大的前三个特征。随后，我们将这些特征映射到它们对应的k-mer区域（图6），详细坐标在表8中提供。我们的分析显示，具有最大预测贡献的区域在CGR图的上右象限有明显的聚类。值得注意的是，在某些属中，贡献最大的位置并不与高k-mer丰度相关，而是与相对低丰度的区域相关，例如葡萄球菌（Staphylococcus）和假单胞菌（Pseudomonas）就是这样的例子。这种现象的发生是因为低频区域表现出跨属的差异性，在统计上更具区分度，而高频区域在多个属中是保守的，因此提供的分类分辨率有限。这些发现强调了PhageCGRNet通过整合k-mer频率和CGR位置信息，有效地捕捉到了低频但高特异性的序列特征，从而提高了噬菌体分类的区分能力。为了验证所识别的低频、高特异性k-mer是否对噬菌体基因组数据库的组成具有鲁棒性，我们进行了重复的子采样实验（n=20），发现贡献最大的k-mer在各个子集中保持一致（重叠率>85%），这证实了检测到的特征并没有受到所使用特定数据集的显著影响。下载：PPT、PowerPoint幻灯片、PNG、更大图像、TIFF、原始图像。表8. 基于SHAP值分析的每个噬菌体属中贡献最大的三个k-mer区域。https://doi.org/10.1371/journal.pcbi.1014130.t008

然后，我们进一步分析了通过SHAP识别出的高贡献序列模式。具体来说，我们检查了噬菌体基因组中最具贡献的序列模式的分布。我们的发现表明，这些模式在对于噬菌体-宿主相互作用至关重要的功能基因中富集，包括尾纤维蛋白基因（介导宿主受体识别）、噬菌体溶菌酶基因（参与细胞壁降解）、尾带测量蛋白基因以及转座酶基因（调节基因组整合），所有这些在先前的研究[41–44]中都有详细记录。PhageCGRNet通过CGR捕捉全局和局部序列模式，同时整合k-mer频率数据，这种双重表示方式融合了多维信息以识别复杂模式。此外，我们从数据集2中选择了在物种水平上有超过100个序列的噬菌体：来自铜绿假单胞菌（Pseudomonas aeruginosa）的117个序列，来自大肠杆菌（Escherichia coli）的226个序列，以及来自粘质芽孢杆菌（Mycolicibacterium smegmatis）的312个序列，总共655个序列。根据第2.4节中的距离计算方法，我们计算了这665个序列的距离矩阵，并使用Neighbor-Joining [45]方法在Mega X中构建了系统发育树，一棵树仅基于k-mer频率（图7(A)），另一棵树结合了频率和位置信息（图7(B)）。比较图7(A)和7(B)可以观察到序列分类的差异。两幅图都将六个铜绿假单胞菌序列错误地分类为大肠杆菌。此外，图7(A)显示几个铜绿假单胞菌序列与粘质芽孢杆菌序列交织在一起。然而，这个问题在图7(B)中不存在，表明使用频率和位置信息作为特征时准确性得到了提高。下载：PPT、PowerPoint幻灯片、PNG、更大图像、TIFF、原始图像。图7. 物种的系统发育树。(A)显示了使用k-mer频率作为特征的构建的系统发育树。(B)显示了使用k-mer频率和位置信息作为特征的构建的系统发育树。https://doi.org/10.1371/journal.pcbi.1014130.g007

对于这六个被错误分类的序列，我们进行了更详细的分析。M11912和M19377是丝状噬菌体，其基因组长度为5,833 bp，明显短于该数据集中其余铜绿假单胞菌噬菌体的平均基因组长度64,868 bp。这种大小差异可能扰乱了k-mer分布模式，并影响了从CGR位置提取特征的效果。此外，三个贡献最大的区域是相同的，这使得特征提取无效，导致了分类错误。另外，序列NC_031274错误地与大肠杆菌的序列JX415536、KF981730、EU330206、NC_027395、JX867715和NC_028903形成了一个簇。根据第4.4节中的方程(13)，我们计算了序列NC_031274与大肠杆菌序列组之间的平均相似度为0.5323。然而，NC_031274与其他铜绿假单胞菌序列之间的平均余弦相似度低于0.5，这导致了分类错误。此外，三个序列HG962376、JQ307387和NC_028879形成了独立的亚组。我们计算了这些序列与其他铜绿假单胞菌序列之间的成对相似度，发现所有值都低于0.3。

2.9. 案例研究 我们选择了Ferry等人[46]报告的临床案例中的三种噬菌体——vB_PaeP_4029、vB_PaeP_4032和vB_PaeP_4034，这些噬菌体成功用于治疗由多重耐药铜绿假单胞菌引起的脊柱感染。这些噬菌体被证实能够特异性地感染并裂解患者的铜绿假单胞菌菌株，包括其持续存在的小菌落变体。这些噬菌体的基因组高度相似，vB_PaeP_4032和vB_PaeP_4034仅相差五个单核苷酸多态性。所有三种噬菌体都属于Limuavirus属，这对模型区分细微序列变异和识别特定基因组模式提出了挑战。我们使用PhageCGRNet根据这三种噬菌体的基因组序列进行宿主预测，并将结果与它们的已知宿主进行比较。结果总结在表9中。结果表明，PhageCGRNet正确预测了所有三种噬菌体的宿主为铜绿假单胞菌，预测概率超过0.9999。这证明了该模型在这种情况下具有高置信度和准确的预测能力。相比之下，RaFAH [11]、vHULK [12]、CL4PHI [25]和PHPGCA [24]也正确预测了宿主，尽管预测概率低于PhageCGRNet。然而，DeepHos t[22]将噬菌体vB_PaeP_4032错误地分类为Gordonia terrae，CHERRY [23]错误地将噬菌体vB_PaeP_4032分类为Pseudomonas syringae，PB-LKS [27]将噬菌体vB_PaeP_4034错误地分配给Synechococcus sp. WH 7805。值得注意的是，CL4PHI [25]报告的是距离值；其他方法提供的是概率值。下载：PPT、PowerPoint幻灯片、PNG、更大图像、TIFF、原始图像。表9. PhageCGRNet对三种噬菌体基因组的预测。https://doi.org/10.1371/journal.pcbi.1014130.t009

此外，我们分析了这三种噬菌体的特征贡献。通过计算SHAP值，我们确定了对预测结果贡献最大的序列模式对应于k-mer GTGCAGG。进一步分析显示，GTGCAGG在编码假定尾纤维蛋白的基因区域内反复出现。我们模型之所以能够准确分类宿主，正是由于它将k-mer频率与CGR位置信息相结合。

3. 讨论与结论 在这项研究中，我们提出了一种名为PhageCGRNet的新深度学习方法，用于预测噬菌体宿主分类。由于噬菌体基因组序列的长度各不相同，我们的方法将这些序列表示为统一大小的三维矩阵。三维矩阵包含两种类型的信息：从序列中获得的k-mer的频率信息以及从CGR图中获得的k-mer的位置信息。首先，我们提取了噬菌体基因组序列的k-mer频率信息。接下来，我们将序列映射到CGR图上。对于每个k-mer，我们找到其在CGR图中结束的子序列的平均位置，作为该k-mer的位置信息。通过这种方式，我们获得了一种新的k-mer度量方法，使得序列可以用这种新的k-mer度量方法来表示。最后，我们将3D特征矩阵作为输入传递给CNN模型，该模型可以处理多个输入通道。CNN使用卷积核扫描并分别求和每个通道的数据，从而结合所有通道的信息进行最终预测。实验结果表明，与其他最先进的方法相比，我们的方法在两个数据集上表现更优。同时，我们在消融实验中发现，CGR提供的位置信息比频率信息更具准确性。CGR衍生的位置信息的有效性可以从其解决生物学上重要序列模式的能力中理解。此外，CGR能够识别生物学上相关的序列模式。例如，Deschavanne等人[39]通过CGR图像观察到四个七字母寡核苷酸（GGCGATC、GCGATCG、CGATCGC和GATCGCC）在Synechocystis sp.中表现出异常高的频率。Robinson等人[47]进一步确认HIP1（5’-GCGATCGC-3’）在蓝细菌Synechococcus中高度富集。这个八聚体可能通过调节mRNA稳定性或翻译终止效率间接调控金属硫蛋白的表达。此外，Scholz等人[48]在Synechocystis sp. (Syn6803)中识别出Cas6作为HIP1，证明了它能够识别并切割pre-crRNA中的回文重复结构以生成成熟的crRNA。这一过程是CRISPR-Cas系统中的关键步骤，因为crRNA引导效应蛋白如Cas9切割目标DNA序列（例如，来自入侵噬菌体或质粒的序列）。这个案例不仅展示了CGR捕捉功能序列模式的能力，还表明这些模式可能具有生物学意义。我们的消融实验表明，结合这两种特征比单独使用任何一种特征都能获得更高的分类准确性。我们结合了这两种特征，获得了比其他先进方法更好的准确性。实际上，捕捉序列相似性是噬菌体-宿主预测的关键。我们结合k-mer频率和位置特征不仅增强了相似性的表示，还编码了具有潜在生物学相关性的基因组模式，这一点得到了宿主相互作用基因中的功能富集的支持（第2.8节）。在未来的工作中，我们计划探索更先进的深度学习架构，如基于Transformer的模型和图神经网络（GNNs），以更好地捕捉噬菌体基因组中的长距离依赖性和结构关系。此外，我们还将研究集成多个单一模型的混合模型（例如，将CNN与循环或注意力机制结合），以增强特征表示和泛化能力。预测的可解释性也将是一个重点，使用注意力可视化等技术来识别生物学上有意义的k-mer模式。最后，我们计划扩展我们的方法以预测更广泛的宿主范围，并整合多模态数据（例如，蛋白质序列或蛋白质结构），以进一步提高分类准确性和实际应用价值。

4. 方法 我们提出的方法包括四个核心组成部分：(i) 输入，(ii) 特征提取，(iii) CNN模块，和(iv) 输出，如图8所示。在输入阶段，我们输入噬菌体的基因组序列。这些序列由四种核苷酸碱基A、T、G和C组成。由于这些序列的长度各不相同，因此从它们中提取特征是必要的。特征提取有助于表示序列，并实现特征矩阵或特征向量的大小统一性。在这里，我们将从序列中获得的k-mer频率信息与从CGR图中获得的k-mer位置信息结合起来。这个过程构成了我们的特征提取阶段。然后将整合的特征矩阵输入CNN模型，以完成噬菌体宿主分类的任务。最后，模型输出分类预测结果。我们在后续章节中提供了每个部分的详细分解。第4.1节涵盖了输入处理和特征提取，第4.2节详细介绍了CNN模型，第4.3节描述了模型训练和输出结果。下载：PPT、PowerPoint幻灯片、PNG、更大图像、TIFF、原始图像。图8. 我们提出的PhageCGRNet的整体过程示意图。首先，每个噬菌体的基因组序列被表示为一个三维矩阵，然后将其输入CNN以进行宿主预测。https://doi.org/10.1371/journal.pcbi.1014130.g008

4.1. 特征矩阵的生成 首先，我们输入噬菌体基因组序列。假设是一个长度为L的DNA序列，其中对于任何，是一个k-mer。首先，我们计算该序列中所有k-mer出现的次数，表示为。然后，根据图9(A)中所示的CGR排列，我们将k-mer的频率记录在一个矩阵中。这个矩阵表示为特征矩阵，其公式如下：下载：PPT、PowerPoint幻灯片、PNG、更大图像、TIFF、原始图像。图9. 特征矩阵构建过程。(A)显示了DNA序列被分解为k-mer（k=2）以生成频率矩阵。（B）显示了序列被投影到CGR图上（红色点代表CGR中每个核苷酸的位置），以生成两个位置矩阵（F2和F3，分别对应x和y坐标）。https://doi.org/10.1371/journal.pcbi.1014130.g009(5)接下来，使用以下迭代函数系统（方程(6)）将DNA序列中的每个核苷酸一一映射到单位平面上，从而生成一个CGR图。(6)其中，我们知道对于DNA序列中的每个k-mer，我们可以识别所有以该k-mer结尾的子序列。随后，我们可以获得这些子序列的CGR。具体来说，我们将以给定k-mer结尾的子序列在CGR图中的位置表示为。需要注意的是，在一个序列中，一个k-mer可能会出现多次。因此，以该k-mer结尾的子序列在CGR图中有多个对应的坐标。为了处理这种情况，我们计算了这些多个坐标点的平均位置 [49]。(7)(8)其中 k-mer 在该序列中出现的次数用 表示。之后，我们分别将 和 视为第二和第三特征，得到矩阵 和 。具体表示如下：(9)(10)为了说明如何生成 3D 矩阵，我们在图 9 中提供了一个示例。假设我们有序列 S = AGTCGTTACA，并且 k-mer 的长度设置为 2。图 9(A) 展示了计算 k-mer 频率的过程。我们将频率信息表示为一个矩阵，其中元素坐标根据 CGR（Color Gradient Representation）进行排列。在图 9(B) 中，我们可以看到将 DNA 序列映射到 CGR 图上的过程。映射完成后，我们计算了相应 k-mer 在 CGR 图上的平均位置，得到 和 的值。由于矩阵 表示 k-mer 的频率，而矩阵 和 是从 CGR 图中点的坐标获得的，因此需要统一度量尺度。为了实现这一点，我们对这三个矩阵进行了归一化。归一化过程如下：(11)此外，我们提供了矩阵 、 和 的详细计算过程。如果序列 S = AGTCGTTACA，k = 2，则序列 S 在 CGR 上的位置坐标如下：A: (0.25, 0.25)，AG: (0.625, 0.625)；AGT: (0.8125, 0.3125)；AGTC: (0.40625, 0.65625)；AGTCG: (0.703125, 0.828125)；AGTCGT: (0.8515625, 0.4140625)；AGTCGTT: (0.92578125, 0.202703125)；AGTCGTTA: (0.462890625, 0.103515625)；AGTCGTTAC: (0.2314453125, 0.5517578125)，AGTCGTTACA: (0.11572265625, 0.27587890625)。因此，平均位置可以计算为：；因此 和 可以计算为：4.2. 卷积神经网络 在构建了噬菌体序列的三维矩阵表示之后，我们使用 CNN 模型来进一步增强特征并分类宿主。具体来说，我们在归一化后将 3D 矩阵输入到 CNN 中。我们设计了具有两个卷积层的 CNN 架构。每个卷积层后面都跟着一个池化层，其中我们采用了 MaxPool 方法 [50]。此外，应用了修正线性单元（ReLU）函数 [51]，定义为 f(x) = max{0, x}，以增强模型在每一层的非线性表达能力。接下来的步骤是这些特征的展平，以便将它们转换为一维向量。之后，我们引入了一个全连接层。这一层将前一层的每个神经元与当前层的所有神经元连接起来，使模型能够学习特征之间的复杂关系。最后，全连接层的输出被送入预测层。在这里，我们使用 Softmax 激活函数 [52] 来计算每个宿主分类的得分。Softmax 函数用于将原始输出值转换为概率分布。然后选择概率得分最高的分类作为模型的预测结果。4.3. 模型训练 在训练过程中，我们使用交叉熵损失函数 [53] 来量化预测的概率分布与宿主分类的真实标签之间的差异。交叉熵损失是分类任务中广泛使用的度量标准。对于多类分类，我们使用分类交叉熵损失，其公式为：(12)其中，N 是类别的数量， 是第 i 类的真实标签， 是第 i 类的预测概率。为了优化我们的 CNN 模型的参数，我们使用了 Adam 优化算法 [54]。我们将该算法的学习率设置为 0.001。最后，模型输出分类预测结果，根据学习到的特征预测每个宿主的类别。4.4. 距离计算 我们将包含 k-mer 频率信息和 CGR 位置信息的 3D 矩阵转置为一个向量。该向量的元素按照特定顺序排列，即 k-mer 频率、CGR 的 x 坐标的平均值以及 CGR 的 y 坐标的平均值。向量转换完成后，我们使用公式 (13) 计算序列相似性，并通过公式 (14) 派生出序列距离。最后，使用 Neighbor-Joining 算法 [45] 构建系统发育树。(13)(14)致谢 本工作得到了中国国家自然科学基金（项目编号 12371088，资助人 ZG Y）和湖南省自然科学基金创新研究群体项目（项目编号 2024JJ1008，资助人 ZG Y）的支持。资助者在研究设计、数据收集与分析、发表决定或手稿准备方面没有发挥作用。参考文献