一种基于机器学习的新冠病毒亚群识别方法与流程

1.本发明属于生物信息领域，具体涉及一种基于机器学习的新冠病毒亚群识别方法。


背景技术：

2.现存的新冠病毒亚群识别主要是在利用比对软件获得病毒突变位点的基础上，根据共同突变将病毒分成不同的亚群。分群的方法主要分为四种，一是通过人工检测的方法把出现相同的高频突变的病毒序列归为一群。这种分析方法有如下两个缺点：1.分群粒度较粗。通过人工方法识别的模式序列，粒度比较粗，即没有找到模式序列的全长，这可能是因为模式序列中存在一些突变频率不那么高的突变，在人工识别的过程中被忽略掉。
3.2.忽略一些潜在的发展中的亚群。由(1)引申出的另一个问题是，病毒群中的亚群，即与病毒群具有相同的主要突变模式，同时又拥有新型的突变的病毒，因为频率不高没有被发现，但是具有形成群体的趋势。这些亚群往往具有时空关系，对于疫情监控，病毒溯源具有重要意义。
4.二是利用进化分析软件根据比对结果构建进化树，根据进化树判断病毒属于哪个种系。这种方法的缺点是当数据量比较多的时候，需要大量的计算资源，普通的计算机没法支撑这样的分析，同时分析也需要大量的时间。其次是当数据量大的时候，进化树的结构会变得相当复杂，难以进行解读。
5.三是利用单体型分析软件根据比对结果构建单体型树。这种方法与第二种类似，同样当数据量比较多的时候，分析和结果解读都会比较困难。而且一种病毒可能会与不同类型的单体型都有关系，因此其具体属于哪种单体型依赖于个体的判断。
6.四是基于一套人为制定的命名规则对病毒进行分群，但是这种方法并不稳定，会间歇出现把具有相同突变模式的病毒划分到不同的群，或者把不同突变模式的病毒划分到同一群中的情况。


技术实现要素：

7.本发明的目的在于提供一种基于机器学习的新冠病毒亚群识别方法，以解决上述背景技术中所提出的问题。
8.为实现上述目的，本发明提供以下技术方案：一种基于机器学习的新冠病毒亚群识别方法，其结构要点在于：包括收集需要分析的病毒基因序列，可以从公共数据库上可以获得大量的新冠病毒序列，然后，根据规则去除质量较低的序列，清理数据后，使用序列比对软件把序列与通用参考序列进行比对，从比对结果可以知道突变出现的位置，接下来把每条病毒序列对应位置上的碱基提取出来进行热编码，从而把每个碱基转变为一个五维的向量，接下来，使用t-sne把矩阵投影到二维平面上，其中t-sne的困惑度为30，然后使用dbscan进行聚类，其中两点之间的最大距离设置为5，最小样本数为50，最后排除掉无法聚
类的噪点，就能得到相互区分开的病毒群，每个群都有特定的模式序列，其中模式序列是通过把突变位点对应回每个病毒群中获得。
9.作为优选的，规则包括：1)未知碱基小于1％，具有长度小于10bp的插入缺失片段，并且唯一突变密度小于等于0.05％；2)未知碱基小于1％，具有长度小于10bp，但是导致编码区移码突变的插入缺失；3)未知碱基小于1％，具有长度大于10bp的插入缺失片段，并且/或突变密度大于0.05％；4)未知碱基比例大于等于1％小于5％；5)长片段的未知碱基，突变密度大于0.05％；6)未知碱基大于等于5％。
10.作为优选的，使用热编码对碱基进行编码，所述的五维的向量中的
“‑”
和“n”表示同一个核苷酸。
11.作为优选的，t-sne为一种基于t分布和条件概率的非线性降维算法。
12.作为优选的，dbscan为一种聚类算法，具有噪声的基于密度的聚类方法。
13.与现有技术相比，本发明具有以下优点1.实现对新冠病毒分群的自动化。传统的进化树和单体型树在得出病毒序列的演化关系后，还需要人工判断哪些病毒属于一个亚群。这里描述的方法结合t-sne和dbscan，实现将突变序列类似的病毒自动聚类，减少人工分群带来的偏差。
14.2.分群粒度更加精细。传统分群方法的粒度较粗，而这种方法能够把病毒亚群分的更细，对疫情的追踪和检测具有重要意义。例如与类似的基于进化树的研究相比，在对gisaid数据库2020年4月之前收集的病毒序列进行分析，我们找到了18个病毒亚群，而对比的研究只找到10个亚群。
15.3.减少计算资源的消耗和缩短分析时间。因为这种方法不需要建立病毒详细的系统发育关系，大大缩减了计算资源和需要分析的时间。在 12cpu，32gb内存的相同的服务器下测试同样的数据集(约3500条病毒数据)，使用这种方法只需要消耗1cpu，1gb内存的资源，在一小时内完成分析任务，而使用传统的进化树方法，在占满所有cpu，和消耗2gb内存的情况下，需要大概5天的时间完成分析。
16.4.快速的病毒溯源。这种方法虽然不能建立系统发育树，但是仍然能够做简单的快速的溯源。例如如果要推断某个群落最早是从哪里爆发的，只需观察该病毒群落中采集时间最早的病毒的采集地点即可。此外，还可以根据病毒群落的模式序列判断群落的演化关系。例如群落a的模式序列为“(p1)a-(p2)g-(p3)t-(p4)t”，而群落b的模式序列为“(p1)a-(p2)g-(p3)t-(p4)t-(p5)c”(p表示突变位点)，因为群落b的模式序列包含了群落a，则可以推断群落b可能是从群落a演化而来。
附图说明
17.图1为本发明对gisaid上的病毒序列进行分群的结果示意图；图2为本发明病毒群的突变信息统计示例；图3为本发明示例中各个病毒群的模式序列示意图。
具体实施方式
18.为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
19.请参阅图1-3，本发明提供一种技术方案，一种基于机器学习的新冠病毒亚群识别方法，包括收集需要分析的病毒基因序列，从gisaid数据库上可以获得大量的新冠病毒序列，然后，根据规则去除质量较低的序列，清理数据后，使用mafftv7.450把序列与通用参考序列nc_045512.2， epi_isl_402125进行比对，以确定突变出现的位置和定位编码区域，蛋白质的注释以nc_045512.2为参比，从比对结果可以知道突变出现的位置，接下来把每条病毒序列对应位置上的碱基提取出来进行热编码，从而把每个碱基转变为一个五维的向量，接下来，使用t-sne把矩阵投影到二维平面上，其中t-sne的困惑度为30，然后使用dbscan进行聚类，其中两点之间的最大距离设置为5，最小样本数为50，最后排除掉无法聚类的噪点，就能得到相互区分开的病毒群，每个群都有特定的模式序列，其中模式序列是通过把突变位点对应回每个病毒群中获得。
20.其中，在本实施例中，所述的规则包括：1)未知碱基小于1％，具有长度小于10bp的插入缺失片段，并且唯一突变密度小于等于0.05％；2)未知碱基小于1％，具有长度小于10bp，但是导致编码区移码突变的插入缺失；3)未知碱基小于1％，具有长度大于10bp的插入缺失片段，并且/或突变密度大于0.05％；4)未知碱基比例大于等于1％小于5％；5)长片段的未知碱基，突变密度大于0.05％；6)未知碱基大于等于5％。
21.其中，在本实施例中，所述的碱基按照以下编码规则进行热编码：1)a编码为(1,0,0,0,0)2)t编码为(0,1,0,0,0)3)g编码为(0,0,1,0,0)4)c编码为(0,0,0,1,0)5)
“‑”
(缺失)和“n”(未知碱基)编码为(0,0,0,0,1)。
22.其中，在本实施例中，所述的t-sne能够使得具有相似突变的病毒序列在低维空间投影后互相靠近，使不相似的序列互相分离开来。t-sne是一种非线性降维算法，能够保持高维数据的结构在低维投影尽可能相似，且能够避免同一区域内的投影点过于密集。也就是说t-sne能够使得具有相似突变的病毒序列在低维空间投影后互相靠近，使不相似的序列互相分离开来。与传统的构建进化树方法相比，这种方法只考虑突变的相似性，无须考虑
重构序列的系统发生顺序，因而大大减少了计算量，从而缩短分析时间和减少计算资源的消耗。
23.其中，在本实施例中，所述的dbscan为一种聚类算法，具有噪声的基于密度的聚类方法。dbscan是一种无须规定聚类数目的算法，能够按照每一类中包含个体的数目以及不同类之间的距离自动聚类。传统方法在得出进化树或者单体型树后，还需要人为地去判断哪些序列是属于同一群体。实施例
24.本实例分析的病毒序列来自于gisaid数据库2020年4月8日之前采集的病毒序列。使用的比对软件是mafftv7.450。分析数据使用的编程语言是python，机器学习方法的实现使用sklearn v0.23.1版本。
25.从gisaid数据库收集2020年4月8日之前，一共4926条新冠病毒全基因组序列以及附属信息，包括采集日期，采集地点等。首先用以下规则对数据进行清洗：1.未知碱基小于1％，具有长度小于10bp的插入缺失片段，并且唯一突变密度小于等于0.05％。
26.2.未知碱基小于1％，具有长度小于10bp，但是导致编码区移码突变的插入缺失。
27.3.未知碱基小于1％，具有长度大于10bp的插入缺失片段，并且(或) 突变密度大于0.05％。
28.4.未知碱基比例大于等于1％小于5％。
29.5.长片段的未知碱基，(或)插入缺失信息(如果有)已知，(或) 突变密度大于0.05％。6.未知碱基大于等于5％。
30.清洗数据后得到3490条高质量的基因组序列。使用比对软件 mafftv7.450比对后一共获得2445个snv(单核苷酸变异)。接着从每条序列中提取2445个位点，位点的位置与snv出现的位置对应。利用这些位点组成一个3490*2445的矩阵，每行表示一个样本，每列表示样本在这个位点的碱基信息。
31.接下来按照以下编码规则对碱基信息进行热编码：1.a编码为(1,0,0,0,0)2.t编码为(0,1,0,0,0)3.g编码为(0,0,1,0,0)4.c编码为(0,0,0,1,0)5.
“‑”
(缺失)和“n”(未知碱基)编码为(0,0,0,0,1)从而得到一个3490*12225(2445*5)的矩阵。
32.使用sklearn包中的t-sne方法对矩阵进行降维，主要参数设置如下；困惑度设置为30；学习速率设置为200；迭代次数设置为1000；距离计算方式设置为欧几里得距离；其他参数设置为默认值。
33.通过降维后得到一个3490*2的矩阵。
34.使用dbscan进行聚类。主要参数设置如下：
点与点的最大距离设置为5；聚类的最小样本数设置为50；距离计算方式设置为欧几里得距离；其他参数设置为默认值。
35.经过聚类和去除噪点后，可以得到18个病毒亚群，使用pythonseaborn包绘制散点图，如附图1所示。
36.统计每个病毒群的变异情况，得到各个病毒群的变异模式序列，例如如附图2所示为群10的统计情况，该群的模式序列为“t(1059)-t(3037)
ꢀ‑
t(14408)-g(23403)-t(25563)”。各病毒群的模式序列如附图3所示。根据突变模式序列的复杂程度以为突变的位置，可以推断病毒的演化路径。例如欧洲变种的新冠病毒的特征为23403位点上的腺嘌呤突变为鸟嘌呤，导致氨基酸由天冬氨酸转变为甘氨酸，即d614g。具有这种突变的病毒亚群为 3，9，10，11，13，14，15，16，17，而这些病毒亚群都共享群3的模式序列“t(3037)-t(14408)-t(23403)”，因此群3可能是这些病毒群的祖先群，而群3 中最早收集到的样本时间为2020-01-28。
37.本发明开发一种基于新冠病毒的突变位点进行自动分群的方法，提高病毒群模式序列的分辨率，从而增强发现病毒亚群的能力，增强分群方法的稳健性，尽可能使具有相同模式序列的病毒划分为同一个群体，具有不同模式序列的病毒群区别开来。
38.对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内，不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
39.此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。