噬菌体宿主属预测方法、装置、设备及存储介质与流程

1.本技术涉及生物测序技术领域，特别是涉及一种噬菌体宿主属预测方法、装置、设备及存储介质。


背景技术：

2.噬菌体作为病毒中数量最多、物种多样性最高的一类，在微生物群落中发挥着重要作用。相较于传统基于培养的病毒发现方法，宏基因组学技术进行一次测序就可获得环境中包括噬菌体在内的所有微生物的序列信息，避免了任何和培养相关的偏差，能够准确快速地确定新噬菌体的序列信息，极大地加快了新噬菌体发现的速度。
3.但正是由于未进行培养，缺乏了病毒的宿主信息，宏基因组学技术无法直接确定发现的新噬菌体的宿主。在微生物群落研究中，无论是针对水平基因转移探索而研究温和噬菌体，还是针对噬菌体疗法而研究烈性噬菌体，噬菌体的宿主信息都是十分关键的。因此，鉴定噬菌体宿主一直是微生物研究领域的热点。
4.宏基因组学技术一般采用以illumina为代表的二代测序方法，二代测序方法直接得到的是以reads为单位的短序列片段，现有技术难以将所有以reads为单位的短序列片段拼接成完整的全基因组序列；并且，现有技术无法实现根据噬菌体的短序列片段，准确预测噬菌体宿主。因此，如何根据噬菌体的短序列片段，直接预测噬菌体宿主成为微生物群落研究中亟待解决的问题。


技术实现要素：

5.本技术实施例提供一种噬菌体宿主属预测方法、装置、设备及存储介质，能够根据噬菌体的短序列片段，直接准确预测噬菌体的宿主属。
6.本技术实施例第一方面提供一种噬菌体宿主属预测方法，所述方法包括：
7.提取待检测噬菌体片段的噬菌体特征；其中，所述噬菌体特征包括第一序列特征和蛋白编码特征；
8.将所述噬菌体特征分别与多个不同候选宿主属的第二序列特征进行配对，得到多个序列特征对；其中，所述第二序列特征是所述候选宿主属的双密码子频率和全基因组5
‑
mer频率；
9.将所述多个序列特征对输入第一预设模型，输出所述多个不同候选宿主属各自的第一得分；其中，第一得分用于表征所述第二序列特征与所述噬菌体特征之间的相关性；
10.根据多个第一得分的分布情况，确定用于进行进一步检测的多个目标候选宿主属；
11.将所述待检测噬菌体片段的第一密码子序列依次输入与每个所述目标候选属下多个原核生物对应的多个第二预设模型，得到每个第二预设模型的第二得分；其中，所述第二得分是密码子马尔可夫链模型与所述第一密码子序列的似然度得分；
12.根据每个第二预设模型的第二得分，计算所述多个不同候选宿主属各自对应的第
三得分；其中，所述第三得分是候选宿主属下分值最高的第二预设模型的第二得分；
13.计算所述第一得分和第三得分的加权平均值，将最大加权平均值对应的候选宿主属确定为目标宿主属。
14.可选地，根据多个第一得分的分布情况，确定用于进行进一步检测的目标候选宿主属，包括：
15.在任意候选宿主属的第一得分大于第一预设阈值的情况下，将该候选宿主属确定为目标候选宿主属；
16.在所述多个第一得分中数值最大者小于第一预设阈值并大于第二预设阈值的情况下，将第一得分大于第三预设阈值的候选宿主属确定为目标候选宿主属；
17.在所述多个第一得分中数值最大者小于第二预设阈值的情况下，将所有候选宿主属确定为目标候选宿主属。
18.可选地，提取待检测噬菌体片段的噬菌体特征，包括：
19.对所述待检测噬菌体片段进行反向补序，得到反向补序片段；
20.分别对所述待检测噬菌体片段和所述反向补序片段进行编码，得到碱基序列特征；
21.对所述待检测噬菌体片段的蛋白质编码区域进行注释；
22.对注释后的待检测噬菌体片段添加经过独热编码后的蛋白编码区域，得到蛋白编码特征；
23.将所述到碱基序列特征和所述蛋白编码特征确定为所述噬菌体特征。
24.可选地，所述方法还包括：
25.根据每个原核生物属与现存噬菌体的侵染关系，确定多个目标原核生物属；
26.将侵染目标原核生物属下任意原核生物的噬菌体确定为噬菌体样本，得到多个噬菌体样本；
27.对所述多个噬菌体样本中的每个噬菌体样本的全基因组进行模拟，生成每个噬菌体样本对应的多个人工短重叠群；
28.对所述多个人工短重叠群中的每个人工短重叠群提取序列特征，得到每个噬菌体样本的第一样本特征；
29.对每个目标原核生物属计算双密码子频率和全基因组5
‑
mer频率，得到每个原核生物属的第二样本特征；
30.对多个第一样本特征和多个第二样本特征一一配对，得到多个训练样本对；
31.利用所述训练样本对待训练模型进行多次训练，得到所述第一预设模型。
32.可选地，对多个第一样本特征和多个第二样本特征一一配对，得到多个训练样本对，包括：
33.在任意噬菌体样本与目标原核生物具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的正向样本对；
34.在任意噬菌体样本与目标原核生物不具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的负向样本对；
35.利用所述训练样本对对待训练模型进行多次训练，得到所述第一预设模型，包括：
36.利用所述正向样本对和所述负向样本对对待训练模型进行多次训练，得到所述第一预设模型。
37.可选地，所述方法还包括：
38.提取每个原核生物蛋白质编码区域的碱基排列次序，得到每个原核生物的第二密码子序列；
39.根据所述第二密码子序列，计算得到每个原核生物的密码子马尔可夫链模型的转移概率矩阵；
40.将所述待检测噬菌体片段的第一密码子序列依次输入与每个所述目标候选属下多个原核生物对应的多个第二预设模型，得到每个第二预设模型的第二得分，包括：
41.依次利用所述多个第二预设模型中每个第二预设模型的转移概率矩阵对所述待检测噬菌体片段的第一密码子序列进行计算，得到每个第二预设模型的第二得分；所述第二预设模型是密码子马尔可夫链模型。
42.本技术实施例第二方面提供一种噬菌体宿主属预测装置，所述装置包括：
43.第一提取模块，用于提取待检测噬菌体片段的噬菌体特征；
44.第一配对模块，用于将所述噬菌体特征分别与多个不同候选宿主属的第二序列特征进行配对，得到多个序列特征对；其中，所述第二序列特征是所述候选宿主属的双密码子频率和全基因组5
‑
mer频率；
45.第一输入模块，用于将所述多个序列特征对输入第一预设模型，输出所述多个不同候选宿主属各自的第一得分；其中，第一得分用于表征所述第二序列特征与所述噬菌体特征之间的相关性；
46.第一确定模块，用于根据多个第一得分的分布情况，确定用于进行进一步检测的目标候选宿主属；
47.第二输入模块，用于将所述待检测噬菌体片段的第一密码子序列依次输入与每个所述目标候选属下多个原核生物对应的多个第二预设模型，得到每个第二预设模型的第二得分；其中，所述第二得分是密码子马尔可夫链模型与所述第一密码子序列的似然度得分；
48.得分计算模块，用于根据每个第二预设模型的第二得分，计算所述多个不同候选宿主属各自对应的第三得分；其中，所述第三得分是候选宿主属下分值最高的第二预设模型的第二得分；
49.第二确定模块，用于计算所述第一得分和第三得分的加权平均值，将最大加权平均值对应的候选宿主属确定为目标宿主属。
50.可选地，所述第一确定模块包括：
51.第一确定子模块，用于在任意候选宿主属的第一得分大于第一预设阈值的情况下，将该候选宿主属确定为目标候选宿主属；
52.第二确定子单元，用于在所述多个第一得分中数值最大者小于第一预设阈值并大于第二预设阈值的情况下，将第一得分大于第三预设阈值的候选宿主属确定为目标候选宿主属；
53.第三确定子单元，用于在所述多个第一得分中数值最大者小于第二预设阈值的情况下，将所有候选宿主属确定为目标候选宿主属。
54.可选地，所述第一提取模块包括：
55.补序子模块，用于对所述待检测噬菌体片段进行反向补序，得到反向补序片段；
56.第一编码子模块，用于分别对所述待检测噬菌体片段和所述反向补序片段进行独热编码，得到碱基序列特征；
57.注释子模块，用于对所述待检测噬菌体片段的蛋白质编码区域进行注释；
58.第二编码子模块，用于对注释后的待检测噬菌体片段添加经过独热编码后的蛋白编码区域，得到蛋白编码特征；
59.第三确定子模块，用于将所述到碱基序列特征和所述蛋白编码特征确定为所述噬菌体特征。
60.可选地，所述装置还包括：
61.第三确定模块，用于根据每个原核生物属与现存噬菌体的侵染关系，确定多个目标原核生物属；
62.第四确定模块，用于将侵染目标原核生物属下任意原核生物的噬菌体确定为噬菌体样本，得到多个噬菌体样本；
63.生成模块，用于对所述多个噬菌体样本中的每个噬菌体样本的全基因组进行模拟，生成每个噬菌体样本对应的多个人工短重叠群；
64.第二提取模块，用于对所述多个人工短重叠群中的每个人工短重叠群提取序列特征，得到每个噬菌体样本的第一样本特征；
65.第一计算模块，用于对每个目标原核生物属计算双密码子频率和全基因组5
‑
mer频率，得到每个原核生物的第二样本特征；
66.第二配对模块，用于对多个第一样本特征和多个第二样本特征一一配对，得到多个训练样本对；
67.训练模块，用于利用所述训练样本对待训练模型进行多次训练，得到所述第一预设模型。
68.可选地，所述第二配对模块包括：
69.第一配对子模块，用于在任意噬菌体样本与目标原核生物具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的正向样本对；
70.第二配对子模块，用于在任意噬菌体样本与目标原核生物不具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的负向样本对；
71.所述训练模块包括：训练子模块，用于利用所述正向样本对和所述负向样本对对待训练模型进行多次训练，得到所述第一预设模型。
72.可选地，所述装置还包括：
73.第三提取模块，用于提取每个原核生物蛋白质编码区域的碱基排列次序，得到每个原核生物的第二密码子序列；
74.第二计算模块，用于根据所述第二密码子序列，计算得到每个原核生物的密码子马尔可夫链模型的转移概率矩阵；
75.所述输入子模块，用于依次利用所述多个第二预设模型中每个第二预设模的转移
概率矩阵对所述待检测噬菌体片段的第一密码子序列进行计算，得到每个第二预设模型的第二得分。
76.本技术实施例第三方面提供一种可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时，实现如本技术第一方面所述的方法中的步骤。
77.本技术实施例第四方面提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时，实现本技术第一方面所述的方法的步骤。
78.本技术实施例从噬菌体由于缺乏trna，在与其宿主协同演化的过程中，会适应其宿主的序列特征(密码子使用偏好、k
‑
mer频率以及gc含量等)的特性出发，将噬菌体的宿主鉴定任务转换为噬菌体片段与其候选宿主之间的侵染关系预测，从而实现有效地以噬菌体的短序列片段作为输入，准确预测噬菌体的宿主属的目的。具体地，本技术实施例提取待检测噬菌体片段的第一序列特征，构建待检测噬菌体片段与每个候选宿主属的配对，使得第一序列特征与每个候选宿主属的第一序列特征组成特征对，将特征对输入预先训练完成的第一预设模型，第一预设模型用于计算噬菌体的短序列片段与候选宿主属的序列特征的相关性，因此第一预设模型输出的第一得分，可以准确表示与第一序列特征组成的特征对的第二序列特征和第一序列特征的相关性，以第二序列特征与第一序列特征的相关性作为判断依据，确定候选宿主属与待检测噬菌体片段是否具有侵染关系，从而确定候选宿主属是否是能够作为进一步检测的多个目标候选宿主属。
79.本技术实施例还构建了每个原核生物的密码子马尔可夫链模型，以第一预设模型输出的第一得分为依据，选择进一步验证原核生物与待检测噬菌体片段的传染关系的多个第二预设模型，将待检测噬菌体片段的密码子分别输入多个第二预设模型，根据多个第二预设模型各自对待检测噬菌体片段的密码子序列的似然度得分，得到多个第二预设模型各自对应的原核生物与待检测噬菌体片段的传染关系，以每个候选宿主属中第二得分最高的原核生物的第二得分分值作为候选宿主属的分值，以再一次验证候选宿主属与待检测噬菌体片段的侵染关系，最后计算候选宿主属的第一得分和根据第二得分确定的候选宿主属的分值的加权平均值，以加权平均值最高的候选宿主属作为待检测噬菌体片段的候选宿主属。
附图说明
80.为了更清楚地说明本技术实施例的技术方案，下面将对本技术实施例的描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本技术的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
81.图1是本技术实施例构建的待训练模型的结构示意图；
82.图2是本技术实施例训练第一预设模型的步骤流程图；
83.图3是本技术一种示例对人工短重叠群提取序列特征的示意图；
84.图4是本技术实施例提取第二样本特征的示意图；
85.图5是本技术实施例中训练样本对的示例图；
86.图6是本技术实施例提出的噬菌体宿主属预测方法的步骤流程图；
87.图7本技术一种示例提出的噬菌体宿主属预测方法的流程图；
88.图8是本技术实施例提出的噬菌体宿主属预测装置的功能结构图。
具体实施方式
89.下面将结合本技术实施例中的附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本技术的一部分实施例，而不是全部的实施例。基于本技术中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本技术保护的范围。
90.现有的检测噬菌体宿主的策略包括：
91.1、计算待测噬菌体与现存噬菌体的相关性，根据现存噬菌体的宿主推断待测噬菌体的宿主。该策略方法对现有的数据库的依赖较大，数据库已记载的噬菌体有限，对于部分新发现的噬菌体，无法有效确定其宿主。微生物群落中存在大量未知的噬菌体，相较于当前数据库中已知噬菌体，序列相关性低。
92.2、多个样品中具有相似丰度分布变化趋势的原核生物可能是该噬菌体的宿主。基于丰度分布的策略需要多个样本来计算每种噬菌体与原核生物之间的相关性，并且噬菌体与其宿主之间的种群动态通常表现出非线性关系，进一步增加了宿主的鉴别难度。
93.3、比对原核生物含有的crispr间隔子和待测噬菌体的序列，基于间隔子的策略仅在噬菌体片段足够长以覆盖crispr区域时才有效，并且仅有部分原核生物均具有crispr系统。
94.鉴于上述问题，本技术实施例提出噬菌体宿主属预测方法，从噬菌体由于缺乏trna，在与其宿主协同演化的过程中，会适应期宿主的序列特征(密码子使用偏好、k
‑
mer频率以及gc含量等)的特性出发，将噬菌体的宿主鉴定任务转换为噬菌体片段与其候选宿主之间的侵染关系预测；根据噬菌体的序列特征与原核生物的序列特征的相关性，确定噬菌体片段和原核生物之间是否存在传染性关系，将与噬菌体片段具有传染性的原核生物属作为噬菌体的宿主属。上述方法避免了序列比对，并且上述密码子使用偏好、k
‑
mer频率以及gc含量等序列特征分布在整个噬菌体基因组，以reads为单位的噬菌体的短序列片段也包含上述信息，因此，本技术实施例提出的噬菌体宿主属预测方法能够有效应对宏基因组测序产生的大量噬菌体短序列片段的宿主预测问题。
95.为了更智能地实施申请人提出的噬菌体宿主属预测方法，使得该方法的应用范围更广，申请人首先构建了待训练模型，并基于噬菌体样本和原核生物样本对待训练模型进行训练，得到第一预设模型。
96.图1是本技术实施例构建的待训练模型的结构示意图，如图1所示，待训练模型由深度学习模型构成，具体采用深度神经网络模型(goolglenet)中的inception模块和常规卷积神经网络构建的。inception模块由多路径的卷积神经网络构成，每条路径使用不同大小的卷积核。通过使用inception模块，模型可以在多个尺度上提取序列特征。
97.待训练模型包括：卷积层(conv)、批归一化层(batchnorm)、最大池化层(maxpool)、均值池化层(averagepool)、全连接层(fc)、随机丢弃层(dropout)、连接层(concat)、输出层(softmax)。
98.图2是本技术实施例训练第一预设模型的步骤流程图，如图2所示，步骤如下：
99.步骤s21：根据每个原核生物属与现存噬菌体的侵染关系，确定多个目标原核生物属。
100.本技术一种示例说明了根据每个原核生物属与现存噬菌体的侵染关系，确定多个目标原核生物属可以的具体实施方式：每个原核生物属与现存噬菌体的侵染关系可以是以下信息的列表：能够侵染芽孢杆菌属的现存噬菌体有142个；能够侵染弧菌属的现存噬菌体有108个；能够侵染志贺氏菌属的现存噬菌体有39个。
101.在特定原核生物属能够被侵染时，将该特定原核生物属确定为目标原核生物属。进一步地为了更好地获得训练第一预设模型的训练集和测试集，本技术实施例在侵染特定原核生物属的现存噬菌体的数量大于2时，将该特定原核生物属确定为目标原核生物属。
102.步骤s22：将侵染目标原核生物属下任意原核生物的噬菌体确定为噬菌体样本，得到多个噬菌体样本。
103.在本技术一种示例中，将芽孢杆菌、染弧菌属以及志贺氏菌属构成集合a，保留可以感染集合a中任意原核生物的噬菌体，得到多个噬菌体样本。
104.在采集样本阶段，可以从ncbi数据库下载现存原核生物和具有宿主注释的噬菌体的完整基因组及注释信息。
105.步骤s23：对所述多个噬菌体样本中的每个噬菌体样本的全基因组进行模拟，生成每个噬菌体样本对应的多个人工短重叠群。
106.为了在训练节点更好地模拟宏基因组学技术得到的噬菌体的短序列片段，使待训练模型适应短序列片段的特征提取和计算，本技术实施例基于噬菌体样本的全基因组，模拟生成人工短重叠群。对噬菌体样本的全基因组进行模拟可以采用基因组学和宏基因组学测序模拟器(metasim)。
107.为了保证测试集中所有测试数据符合实际应用中噬菌体的短序列片段第一次被发现的特点，本技术一种具体实施方式按照噬菌体样本首次发布的日期，将噬菌体样本的全基因组划分为训练集和测试集。例如，将2015年之前发布的噬菌体样本的全基因组划分为训练待训练模型的训练集，将2015年之后发布的噬菌体样本的全基因组划分为测试待训练模型的测试集，并调整噬菌体样本的全基因组在训练集和训练集中的分布，以使训练数据测试数据比例合理。示例地，可以侵染芽孢杆菌属的现存噬菌体有142个，2015年之前的噬菌体有21个，此时测试集和训练集中的样本分布不合理，将训练集中可以侵染芽孢杆菌属的数据调整为90个噬菌体样本的全基因组。
108.进一步地为了确保训练集和独立测试集中没有相同或接近相同的噬菌体，本技术在一种具体的实施方式中，在调整噬菌体样本的全基因组在训练集和训练集中的分布时，遵循以下原则：1)相同的噬菌体全基因组不会同时存在于训练集和测试集中；2)同一研究中不同的噬菌体全基因组必须同时存在于训练集或测试集中。
109.为了更好地适应实际应用中噬菌体短序列片段长短可能不一致的情况，本技术一种实施方式分别模拟生成了三组不同长度的人工短重叠群。生成的多个人工短重叠群可以划分为以下三种尺寸：100
‑
400bp，401
‑
800bp以及801
‑
1,200bp。具体地，测试集和训练集中的人工短重叠群都可以划分为100
‑
400bp，401
‑
800bp以及801
‑
1,200bp三种尺寸。相应地，本技术实施方式分别构建了分别对应上述三种尺寸的待训练模型。例如，对应801
‑
1,200bp的待训练模型的输入层与801
‑
1,200bp匹配。
110.步骤s24：对所述多个人工短重叠群中的每个人工短重叠群提取序列特征，得到每个噬菌体样本的第一样本特征。
111.第一样本特征是指噬菌体样本的序列特征。
112.对单个人工短重叠群提取序列特征具体可以是，根据人工短重叠群的碱基信息，对人工短重叠群采用独热(one
‑
hot)编码操作，得到第一矩阵；再对人工短重叠群进行反向补序，对反向补序后的人工短重叠群采用独热编码操作，得到第二矩阵，对第一矩阵和第二矩阵进行拼接，得到第三矩阵。再通过ncbi中噬菌体基因组的注释信息注释人工短重叠群的蛋白质编码区域(cds区域)，对人工短重叠群的cds信息进行独热编码操作，得到第四矩阵。将第三矩阵和第四矩阵作为第一样本特征。
113.图3是本技术一种示例对人工短重叠群提取序列特征的示意图，如图3所示，本技术一种示例中，人工短重叠群的碱基信息是acgctattgcaccg；采用[1,0,0,0]表示碱基a，采用[0,1,0,0]表示碱基c，采用[0,0,1,0]表示碱基g，采用[0,0,0,1]表示碱基t，得到4
×
len_group的矩阵。针对反向补序后的人工短重叠群，采用上述的编码方式，得到另一个4
×
len_group的矩阵。拼接上述两个4
×
len_group的矩阵，得到图3所示的第一样本特征中的输入1，也就是输入1是拼接上述两个4
×
len_group的矩阵，得到的4
×
2len_group的矩阵。针对人工短重叠群中注释的cds区域，采用[1,0]表示非编码区，[0,1]表示编码区，得到图3所示的第一样本特征中的输入2，输入2是2
×
2len_group矩阵。第一样本特征包括输入1和输入2。
[0114]
步骤s25：对每个目标原核生物属计算双密码子频率和全基因组5
‑
mer频率，得到每个原核生物的第二样本特征；
[0115]
图4是本技术实施例提取第二样本特征的示意图，如图4所示，按照原核生物的属对多个目标原核生物进行聚类，计算每个类内原核生物基因组cds区域的双密码子频率和全基因组5
‑
mer频率。特定属的双密码子频率是64
×
64矩阵，作为输入3；特定属的5
‑
mer频率是1024维向量，作为输入4。
[0116]
步骤s26：对多个第一样本特征和多个第二样本特征一一配对，得到多个训练样本对。
[0117]
对多个第一样本特征和多个第二样本特征一一配对具体方法是：在任意噬菌体样本与目标原核生物具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的正向样本对。
[0118]
在任意噬菌体样本与目标原核生物不具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的负向样本对。
[0119]
得到正向样本对和负向样本对后，对正向样本打上标签1，标签1表示第一样本特征提取自可以侵染第二样本特征所属原核生物的噬菌体；对负向样本对打上签0，标签0表示第一样本特征提取自不可以侵染第二样本特征所属原核生物的噬菌体。
[0120]
图5是本技术实施例中训练样本对的示例图，如图5所示，训练样本对包括一个第一样本特征、一个第二样本特征以及第一样本特征和一个第二样本特征的关系。
[0121]
步骤s27：利用所述训练样本对待训练模型进行多次训练，得到所述第一预设模型。具体地，利用所述正向样本对和所述负向样本对对待训练模型进行多次训练，得到所述
第一预设模型。
[0122]
本技术实施例将特定噬菌体的序列特征和特定原核生物的序列特征组成样本对，那么对于任意的样本对，特定噬菌体的序列特征和特定原核生物的序列特征只可能存在相关关系，或者不存在相关关系两种结果，具体在特定噬菌体与特定原核生物存在传染关系时，即特定原核生物是特定噬菌体的宿主时，特定噬菌体的序列特征和特定原核生物的序列特征存在相似关系，在特定噬菌体与特定原核生物不存在传染关系时，特定噬菌体的序列特征和特定原核生物的序列特征存在不相似关系。利用样本对训练待训练模型，使得待训练模型具有能够准确计算新发现的噬菌体的序列特征和特定原核生物的序列特征之间的相关性，以该相关性作为依据，判断该特定原核生物是否是新发现噬菌体的宿主。将判断噬菌体宿主属的任务转换为判断噬菌体与原核生物之间是否存在侵染关系的二分类任务，提高了判断准确率。
[0123]
继续参考图1，具体在训练待训练模型得到第一预设模型的过程中，将输入1和输入2共同输入卷积层1，输入3输入卷积层4，输入4输入全连接层1。待训练模型输出一个打分，根据标签和待训练模型输出的打分之间的损失函数，调整待训练模型的参数，直到待训练模型能准确输出样本对中第一样本特征和第二样本特征的相关性。进一步地，本技术实施例还可以利用验证集中的测试数据对第一预设模型进行测试，进一步优化第一预设模型的参数。
[0124]
本技术上述实施例主要提出了基于样本对训练第一预设模型的训练过程，并最终得到第一预设模型。以下本技术将着重介绍噬菌体宿主属预测方法，并是示意性地介绍如何将第一预设模型用于噬菌体宿主属预测方法中。
[0125]
图6是本技术实施例提出的噬菌体宿主属预测方法的步骤流程图，如图6所示，步骤如下：
[0126]
步骤s61：提取待检测噬菌体片段的噬菌体特征；其中，所述噬菌体特征包括第一序列特征和蛋白编码特征。
[0127]
待检测噬菌体片段是通过宏基因组学技术进行一次测序，得到的环境中所有微生物的序列信息中的特定噬菌体的短序列片段。第一序列特征是短序列片段的碱基序列特征。
[0128]
本技术另一种实施例提出提取第一序列特征的具体方法，方法包括以下步骤：
[0129]
步骤s61
‑
1：对所述待检测噬菌体片段进行反向补序，得到反向补序片段；对待检测噬菌体片段进行反向补序的方法，与训练第一预设模型阶段提取第一样本特征时，对人工短重叠群进行反向补序的方法相同，本技术实施例在此不作赘述。
[0130]
步骤s61
‑
2：分别对所述待检测噬菌体片段和所述反向补序片段进行独热编码，得到碱基序列特征。对待检测噬菌体片段和反向补序片段进行编码的方法，与训练第一预设模型阶段对人工短重叠群采用独热(one
‑
hot)编码操作相同，本技术实施例在此不作赘述。
[0131]
步骤s61
‑
3：对所述待检测噬菌体片段的蛋白质编码区域进行注释。对待检测噬菌体片段的蛋白质编码区域(cds区域)进行注释的方法，与训练第一预设模型阶段对人工短重叠群进行注释的方法不同的是，本技术实施例在具体应用第一预设模型检测待检测噬菌体片段的过程中，待检测噬菌体片段没有相关注释信息，因此在本技术一种示例中采用基因预测工具prodiga对待检测噬菌体片段注释cds区域。
[0132]
步骤s61
‑
4：对注释后的待检测噬菌体片段进行编码，得到蛋白编码特征。对注释后的待检测噬菌体片段进行编码的方法，与训练第一预设模型阶段对人工短重叠群的cds信息进行独热编码操作相同，本技术实施例在此不作赘述。
[0133]
步骤s61
‑
5：将所述到碱基序列特征和所述蛋白编码特征确定为所述噬菌体特征。
[0134]
步骤s62：将所述噬菌体特征分别与多个不同候选宿主属的每个第二序列特征进行配对，得到多个序列特征对；其中，所述第二序列特征是所述候选宿主属的双密码子频率和全基因组5
‑
mer频率。
[0135]
候选宿主属是预先准备的原核生物属，示例地，芽孢杆菌、染弧菌属以及志贺氏菌属等都可以作为候选宿主属。提取候选宿主属的双密码子频率和全基因组5
‑
mer频率的方法，与训练第一预设模型阶段提取目标原核生物的第二样本特征的方法相同，本技术实施例在此不作赘述。
[0136]
本技术实施例提取候选宿主属中多个原核生物的双密码子频率和全基因组5
‑
mer频率，作为候选宿主属的序列特征，以检测待检测噬菌体片段中是否出现与原核生物的双密码子频率和全基因组5
‑
mer频率相似的特征，从而计算待检测噬菌体片段与候选宿主属的序列特征的相似性。
[0137]
步骤s63：将所述多个序列特征对输入第一预设模型，输出所述多个不同候选宿主属各自的第一得分；其中，第一得分用于表征所述第二序列特征与所述噬菌体特征的之间的相关性。
[0138]
在本技术一种示例中，可以根据第一得分，从多个不同候选宿主属中确定所述待检测噬菌体片段的目标宿主属。具体方法如下：将从待检测噬菌体片段atcctggtgcaccg提取的第一序列特征1，分别与芽孢杆菌的第二序列特征1、染弧菌属的第二序列特征2以及志贺氏菌属的第二序列特征3构成序列特征对。序列特征对1是：(第一序列特征1，第二序列特征1)，序列特征对2是：(第一序列特征1，第二序列特征2)，序列特征对3是：(第一序列特征1，第二序列特征3)。将序列特征对1、序列特征对2以及序列特征对3输入第一预设模型，第一预设模型分别针对序列特征对1、序列特征对2以及序列特征对3输出得分：序列特征对1
‑
0.8、序列特征对2
‑
0.85以及序列特征对3
‑
0.2。由此可见，待检测噬菌体片段atcctggtgcaccg与芽孢杆菌序列特征的相关性是0.8，待检测噬菌体片段atcctggtgcaccg与染弧菌属序列特征的相关性是0.85，待检测噬菌体片段atcctggtgcaccg与志贺氏菌属序列特征的相关性是0.2，从而确定染弧菌属是待检测噬菌体片段的目标宿主属。
[0139]
由于目前数据库中具有宿主注释的噬菌体相对有限，并且许多具有宿主注释的噬菌体的宿主分布在一小部分原核生物属的上。由于机器学习算法在很大程度上依赖于训练数据集，因此任何机器学习算法都难以处理没有任何相关知识的全新数据。因此，为了更加准确地预测噬菌体宿主属，本技术另一种实施例根据第一预设模型输出的候选宿主属的得分，选择可能是目标宿主属下的目标原核生物，以马尔可夫链模型的似然度得分作为依据，判断目标原核生物的马尔可夫链模型进一步验证待检测噬菌体片段是否能够适应目标原核生物的马尔可夫链模型，即验证待检测噬菌体片段的序列特征的状态转换与目标原核生物的序列特征的状态转换是否高度相似。
[0140]
本技术实施例首先针对每个现存的原核生物，训练了密码子马尔可夫链模型，即第二预设模型(hophage
‑
s)。
[0141]
首先，提取每个原核生物蛋白质编码区域的碱基排列次序，得到每个原核生物的密码子序列。示例地，碱基序列atgggctactaa的密码子序列为atg、ggc、tac、taa。
[0142]
再根据所述密码子序列，计算得到每个原核生物的密码子马尔可夫链模型的转移概率矩阵；
[0143]
对于候选宿主属下的每个原核生物，根据dna序列数据库(genebank)中的注释信息提取每个原核生物全基因组的cds区域。以表示原核生物i上的第n个cds区域，n
i
表示原核生物i中cds的总数，x1x2......x
k
表示一段特定密码子序列，以#x1x2......x
k
表示原核生物i上的第n个cds区域中特定密码子序列的数量；原核生物i的密码子马尔可夫链模型的转移概率矩阵p
i
可以通过(1)式计算得到：
[0144][0145]
基于第二预设模型的转移概率矩阵p
i
计算待检测噬菌体片段与密码子马尔可夫链模型的似然度得分，以检验待检测噬菌体片段中密码子的状态转移特征是否能够匹配原核生物的密码子马尔可夫链模型的转移概率矩阵。
[0146]
本技术实施例提出的噬菌体宿主属预测方法还包括步骤：
[0147]
步骤s64：根据多个第一得分的分布情况，确定用于进行进一步检测的目标候选宿主属。
[0148]
根据多个第一得分的分布情况确定用于进行检测的目标得分对应的目标候选宿主属包括以下具体方法：在任意候选宿主属的第一得分大于第一预设阈值的情况下，将该候选宿主属确定为目标候选宿主属；在所述多个第一得分中数值最大者小于第一预设阈值并大于第二预设阈值的情况下，将第一得分大于第三预设阈值的候选宿主属确定为目标候选宿主属；在所述多个第一得分中数值最大者小于第二预设阈值的情况下，将所有候选宿主属确定为目标候选宿主属。
[0149]
步骤s65：将所述待检测噬菌体片段的第一密码子序列依次输入与每个所述目标候选属下多个宿主生物对应的多个第二预设模型，得到每个第二预设模型的第二得分；其中，所述第二得分是密码子马尔可夫链模型与所述第一密码子序列的似然度得分。
[0150]
步骤s66：根据每个第二预设模型的第二得分，计算所述多个不同候选宿主属各自对应的第三得分；其中，所述第三得分是候选宿主属下分值最高的第二预设模型的第二得分；
[0151]
步骤s67：计算所述第一得分和第三得分的加权平均值，将最大加权平均值对应的候选宿主属确定为目标宿主属。
[0152]
将所述待检测噬菌体片段的密码子序列依次输入与每个所述目标候选属下多个原核生物对应的多个第二预设模型，得到每个第二预设模型的第二得分，包括：
[0153]
依次利用所述多个第二预设模型中每个第二预设模的转移概率矩阵对所述待检测噬菌体片段的密码子序列进行计算，得到每个第二预设模型的第二得分。
[0154]
在hophage
‑
s中设置k＝2，得到密码子马尔可夫链模型的转移概率是4096
×
64的矩阵。采用(2)式计算第二得分score
i
：
[0155][0156]
其中，表示待检测噬菌体片段第m个cds区域的第j个密码子，m表示待检测噬菌体片段上cds区域的总个数，l
m
表示第m个cds区域上的密码子数目，p
i
是转移概率矩阵。
[0157]
在待检测噬菌体片段没有cds区域时，将直接从待检测噬菌体片段的六个相位中提取六个密码子序列，作为待检测噬菌体片段的密码子特征。
[0158]
图7本技术一种示例提出的噬菌体宿主属预测方法的流程图，如图7所示，对于获取的待检测噬菌体片段，先使用基因预测工具prodigal对其注释cds区域，然后根据注释到的cds区域，提取待检测噬菌体片段的噬菌体特征，构建噬菌体特征和候选宿主属的第二序列特征的特征对。第一序列特征和第二序列特征输入第一预设模型hophage
‑
g，输出第一得分。图7中max(score)为最大第一得分。
[0159]
根据第一得分，选择密码子马尔科夫链模型，在第一得分大于0.8时，hophage
‑
g的结果非常可靠，在hophage
‑
s中仅保留hophage
‑
s里来自hophage
‑
g内打分大于0.8的属的候选宿主属中的多个密码子马尔科夫链模型。候选宿主属中的多个密码子马尔科夫链模型是：候选宿主属下每个原核生物的密码子马尔科夫链模型。在第一得分位于0.40.8的区间时，hophage
‑
g结果的可靠性较高，在hophage
‑
s中，保留hophage
‑
s里来自hophage
‑
g内打分大于0.25的候选宿主属的多个密码子马尔科夫链模型。在第一得分小于0.4时，hophage
‑
g的深度学习模型打分的可靠性很低，将所有候选宿主属下每个原核生物的密码子马尔科夫链模型保留，作为进一步对检测噬菌体片段进行计算的第二预设模型。
[0160]
将待检测噬菌体片段输入确定的多个第二预设模型，即多个原核生物的密码子马尔科夫链模型，确定最高似然度得分对应的密码子马尔科夫链模型1，以密码子马尔科夫链模型1作为，该密码子马尔科夫链模型对应的原核生物所在属的第三得分，计算所述第一得分和第三得分的加权平均值，将最大加权平均值对应的候选宿主属确定为目标宿主属。
[0161]
还可以进一步将hophage
‑
s获得的所有分数均归一化到[0，score_gmax]。因此，hophage
‑
s中一个原核生物属内最高得分和hophage
‑
g中该属的得分进行的加权平均值将作为该属的综合得分。
[0162]
本技术实施例从噬菌体由于缺乏trna，在与其宿主协同演化的过程中，会适应期宿主的序列特征(密码子使用偏好、k
‑
mer频率以及gc含量等)的特性出发，将噬菌体的宿主鉴定任务转换为噬菌体片段与其候选宿主之间的侵染关系预测，从而实现有效地以噬菌体的短序列片段作为依据，预测噬菌体的宿主属的目的。具体地，本技术实施例提取待检测噬菌体片段的第一序列特征，构建待检测噬菌体片段与每个候选宿主属的配对，使得第一序列特征与每个候选宿主属的第一序列特征组成特征对，将特征对输入预先训练完成的第一预设模型，第一预设模型用于计算噬菌体的短序列片段与候选宿主属的序列特征的相关性，因此第一预设模型输出的第一得分，可以准确表示与第一序列特征组成特征对的第二序列特征和第一序列特征的相关性，以第二序列特征与第一序列特征的相关性作为判断依据，确定候选宿主属与待检测噬菌体片段是否具有侵染关系，从而确定候选宿主属是否是待检测噬菌体片段的目标宿主属。
[0163]
同时，本技术实施例还构建了每个原核生物的密码子马尔可夫链模型，以第一预
设模型输出的第一得分为依据，选择进一步验证原核生物与待检测噬菌体片段的传染关系的多个第二预设模型，将待检测噬菌体片段的密码子分别输入多个第二预设模型，根据多个第二预设模型各自对待检测噬菌体片段的密码子序列的似然度得分，得到多个第二预设模型各自对应的原核生物与待检测噬菌体片段的传染关系，以最大似然度得分对应的原核生物所在的目标候选属确定为目标宿主属，进一步验证了待检测噬菌体片段的目标宿主属。
[0164]
图7中输出综合打分最高的宿主属的依据是(3)式，其中，score_genus
n
是特定宿主属的综合打分，weigt_s是hophage
‑
s的权重，score_g_genus
n
是hophage
‑
g对特定宿主属的打分，score_s in genus
n
是hophage
‑
s对特定宿主属的打分。
[0165]
基于同一发明构思，本技术实施例提供一种噬菌体宿主属预测装置。图8是本技术实施例提出的噬菌体宿主属预测装置的功能结构图。如图8所示，该装置包括：
[0166]
第一提取模块81，用于提取待检测噬菌体片段的噬菌体特征；其中，所述噬菌体特征包括第一序列特征和蛋白编码特征；
[0167]
第一配对模块82，用于将所述噬菌体特征分别与多个不同候选宿主属的第二序列特征进行配对，得到多个序列特征对；其中，所述第二序列特征是所述候选宿主属的双密码子频率和全基因组5
‑
mer频率；
[0168]
第一输入模块83，用于将所述多个序列特征对输入第一预设模型，输出所述多个不同候选宿主属各自的第一得分；其中，第一得分用于表征所述第二序列特征与所述噬菌体特征之间的相关性；
[0169]
第一确定模块84，用于根据多个第一得分的分布情况，确定用于进行进一步检测的目标候选宿主属；
[0170]
第二输入模块85，用于将所述待检测噬菌体片段的第一密码子序列依次输入与每个所述目标候选属下多个原核生物对应的多个第二预设模型，得到每个第二预设模型的第二得分；其中，所述第二得分是密码子马尔可夫链模型与所述第一密码子序列的似然度得分；
[0171]
得分计算模块86，用于根据每个第二预设模型的第二得分，计算所述多个不同候选宿主属各自对应的第三得分；其中，所述第三得分是候选宿主属下分值最高的第二预设模型的第二得分；
[0172]
第二确定模块87，用于计算所述第一得分和第三得分的加权平均值，将最大加权平均值对应的候选宿主属确定为目标宿主属。
[0173]
可选地，所述第一确定模块包括：
[0174]
第一确定子模块，用于在任意候选宿主属的第一得分大于第一预设阈值的情况下，将该候选宿主属确定为目标候选宿主属；
[0175]
第二确定子单元，用于在所述多个第一得分中数值最大者小于第一预设阈值并大于第二预设阈值的情况下，将第一得分大于第三预设阈值的候选宿主属确定为目标候选宿主属；
[0176]
第三确定子单元，用于在所述多个第一得分中数值最大者小于第二预设阈值的情况下，将所有候选宿主属确定为目标候选宿主属。
[0177]
可选地，所述第一确定子模块包括：
[0178]
第一确定子单元，用于在任意候选宿主属的第一得分大于第一预设阈值的情况下，将该候选宿主属确定为目标候选宿主属；
[0179]
第二确定子单元，用于在所述多个第一得分中数值最大者小于第一预设阈值并大于第二预设阈值的情况下，将第一得分大于第三预设阈值的候选宿主属确定为目标候选宿主属；
[0180]
第三确定子单元，用于在所述多个第一得分中数值最大者小于第二预设阈值的情况下，将所有候选宿主属确定为目标候选宿主属。
[0181]
可选地，所述第一提取模块包括：
[0182]
补序子模块，用于对所述待检测噬菌体片段进行反向补序，得到反向补序片段；
[0183]
第一编码子模块，用于分别对所述待检测噬菌体片段和所述反向补序片段进行独热编码，得到碱基序列特征；
[0184]
注释子模块，用于对所述待检测噬菌体片段的蛋白质编码区域进行注释；
[0185]
第二编码子模块，用于对注释后的待检测噬菌体片段添加经过独热编码后的蛋白编码区域，得到蛋白编码特征；
[0186]
第五确定子模块，用于将所述到碱基序列特征和所述蛋白编码特征确定为所述噬菌体特征。
[0187]
可选地，所述装置还包括：
[0188]
第二确定模块，用于根据每个原核生物属与现存噬菌体的侵染关系，确定多个目标原核生物属；
[0189]
第三确定模块，用于将侵染目标原核生物属下任意原核生物的噬菌体确定为噬菌体样本，得到多个噬菌体样本；
[0190]
生成模块，用于对所述多个噬菌体样本中的每个噬菌体样本的全基因组进行模拟，生成每个噬菌体样本对应的多个人工短重叠群；
[0191]
第二提取模块，用于对所述多个人工短重叠群中的每个人工短重叠群提取序列特征，得到每个噬菌体样本的第一样本特征；
[0192]
第一计算模块，用于对每个目标原核生物属计算双密码子频率和全基因组5
‑
mer频率，得到每个原核生物的第二样本特征；
[0193]
第二配对模块，用于对多个第一样本特征和多个第二样本特征一一配对，得到多个训练样本对；
[0194]
训练模块，用于利用所述训练样本对待训练模型进行多次训练，得到所述第一预设模型。
[0195]
可选地，所述第二配对模块包括：
[0196]
第一配对子模块，用于在任意噬菌体样本与目标原核生物具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的正向样本对；
[0197]
第二配对子模块，用于在任意噬菌体样本与目标原核生物不具有侵染关系时，将该噬菌体样本的第一样本特征和所述目标原核生物的第二样本特征进行配对，得到所述目标原核生物属的负向样本对；
[0198]
所述训练模块包括：训练子模块，用于利用所述正向样本对和所述负向样本对对
待训练模型进行多次训练，得到所述第一预设模型。
[0199]
可选地，所述装置还包括：
[0200]
第三提取模块，用于提取每个原核生物蛋白质编码区域的碱基排列次序，得到每个原核生物的第二密码子序列；
[0201]
第二计算模块，用于根据所述第二密码子序列，计算得到每个原核生物的密码子马尔可夫链模型的转移概率矩阵；
[0202]
所述输入子模块，用于依次利用所述多个第二预设模型中每个第二预设模的转移概率矩阵对所述待检测噬菌体片段的第一密码子序列进行计算，得到每个第二预设模型的第二得分。
[0203]
基于同一发明构思，本技术另一实施例提供一种可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如本技术上述任一实施例所述的噬菌体宿主属预测方法中的步骤。
[0204]
基于同一发明构思，本技术另一实施例提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行时实现本技术上述任一实施例所述的噬菌体宿主属预测方法中的步骤。
[0205]
对于装置实施例而言，由于其与方法实施例基本相似，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
[0206]
本说明书中的各个实施例均采用递进或说明的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。
[0207]
本领域内的技术人员应明白，本技术实施例的实施例可提供为方法、装置、或计算机程序产品。因此，本技术实施例可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd
‑
rom、光学存储器等)上实施的计算机程序产品的形式。
[0208]
本技术实施例是参照根据本技术实施例的方法、装置、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理终端设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理终端设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0209]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理终端设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0210]
这些计算机程序指令也可装载到计算机或其他可编程数据处理终端设备上，使得在计算机或其他可编程终端设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程终端设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0211]
尽管已描述了本技术实施例的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本技术实施例范围的所有变更和修改。
[0212]
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。
[0213]
以上对本技术所提供的一种噬菌体宿主属预测方法、装置、设备及存储介质，进行了详细介绍，以上实施例的说明只是用于帮助理解本技术的方法及其核心思想；同时，对于本领域的一般技术人员，依据本技术的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本技术的限制。