使用机器学习对核酸序列进行分类的方法

技术领域

本发明涉及生物技术领域，具体涉及使用机器学习对核酸序列进行分类的方法。

背景技术

目前对核酸序列进行分类多采用序列比对的方法。首先获取带有明确标签的核酸序列构建数据库，然后将待检测的核酸序列与数据库中的核酸序列进行逐一比对，找出与其相似性最高的序列，最后此相似性最高的序列所属标签将赋予待检测序列作为它的类别。

序列比对方法虽然使用十分广泛，但是其局限性也十分明显。由于目前数据库均是基于已知分类序列进行构建，通过序列比对方法对新出现的核酸序列则不能正确地进行分类；同时在基因序列测定过程中，受技术和资金所限，获得的核酸序列可能只是其完整基因组序列的一部分，有时长度甚至短于1000bp，而序列比对方法对于短序列的区分度更低；随着数据库中参考序列数目的不断增加，检测目标序列的运算量也会不断加大，甚至在普通计算机上难以完成。

发明内容

为了解决上述问题，本发明提供一种使用机器学习对核酸序列进行分类的方法，所述方法包括以下步骤：S0：收集带有准确分类标签的核酸序列构建核酸序列数据库；S1：将核酸序列数据库中所有核酸序列分割成等长度的短序列，从而形成短核酸序列均衡样本，该短序列的长度为4bp～核酸序列数据库中最短序列的长度，短核酸序列均衡样本适合于机器学习，短核酸序列均衡样本是指最多分类标签数目的短核酸序列数量N

在一种实施方式中，本发明方法还包括步骤S7:展示待检测核酸序列的分类结果。

在一种实施方式中，S7中展示部分为每条序列的ID和该ID序列可能的分类标签；并显示该分类标签的可靠性分值，分值越高则表明这条ID序列为此分类标签可靠性越高。

在一种实施方式中，所述S1包括以下步骤：

S1-1：判断核酸序列数据库中所有核酸序列是否均衡；

S1-2-1：如果均衡，则使用相同分割步长将核酸序列数据库中所有核酸序列分割成等长度的短序列，所述相同分割步长是分割时后移长度，该后移长度不大于所述短序列的长度，分割核酸序列数据库中所有核酸序列获得短序列样本；

S1-2-2：如果不均衡，则对核酸序列数据库中每类标签核酸序列使用各自分割步长将每类标签核酸序列分割成等长度的短序列，每类标签核酸序列分割步长为分割各类标签核酸序列时后移长度，该后移长度不大于所述短序列的长度，分割核酸序列数据库中所有核酸序列获得短序列样本；

S1-3：判断S1-2-1或者S1-2-2获得短核酸序列的样本是否为均衡并适合于机器学习的样本，如果是，则停止S1-2-1或者S1-2-2，如果不是，则重复S1-2-1或者S1-2-2，直到获得均衡并适合于机器学习的短核酸序列的样本。

在一种实施方式中，所述S1-2-2中，包括以下步骤：(1)获得核酸序列数据库中带有准确分类标签的核酸序列中每类核酸序列数目与特定值的比率，特定值的选择范围在核酸序列数据库中核酸序列中最多标签数目的核酸序列数量N

在一种实施方式中，S2包括：S2-1：使用软件提取核酸序列数据库中所有的Kmer组成集合；和，S2-2：构建Kmer字典，将S2-1中的所有Kmer按字母顺序从小到大排序，并且排序好后的Kmer中的第一位置插入表示字典中不存在的Kmer。

在一种实施方式中，所述字典的查询方式为Kmer:Kmer所在字典中的位置。

在一种实施方式中，S2-1中软件是jellyfish软件。

在一种实施方式中，所述S4包括以下步骤：S4-1：搭建一个含有输入层、Embedding层、最大池化层、平均池化层、中间合并层、全连接层一、全连接层二、输出层的模型；S4-2：把S3得到的数字信号输入到S4-1搭建的模型中进行训练，直到该神经网络成为对输入核酸序列进行多分类的神经网络分类器，分类标签为S0核酸数据库中的分类标签；S4-3：把S4-2中训练好的模型保存，用于后续对待检测序列进行分析。

在一种实施方式中，使用S3中50％-90％数字信号作为训练集，S3中10％-50％数字信号作为测试集，将训练集输入S4-2中搭建的神经网络，并对该网络进行反复训练，并利用测试集不断测试，直到该网络的分类表现达到成为对输入核酸序列进行多分类的神经网络分类器的要求，其分类标签为核酸数据库中核酸序列的分类标签。

在本发明方法中，在步骤S1中，根据数据库中各类别核酸的数目设定相应的步长对各类别的核酸进行分割，从而获得短序列均衡样本，使其适合于机器学习并得到可靠的神经网络分类器，短序列的长度可根据后续待检测样本的长度灵活选择。

在S2步骤中，使用Kmer提取软件直接提取数据库中所有的Kmer组成字典，可以减少后续步骤中与字典对应的Embedding数组的运算量，减少模型的存储空间，增加模型的可移植性。

在S4步骤中，训练一个核酸序列分类器并保存，可以快速将待检测核酸序列与数据库中的核酸序列标签进行对应分类。

在S7步骤中，展示部分为每条序列的ID和该ID序列多种可能的分类标签；并显示可能分类标签的可靠性分值，分值越高则表明这条ID序列为此分类标签可靠性越高，为研究该ID序列与多类标签之间的关系提供参考。

通过本发明的方法，能够准确快速地将待检测核酸序列与数据库中的核酸序列标签进行对应分类，并展示每个标签的可靠性分值，为后续研究提供参考，增加了模型的使用价值。

当数据库中核酸序列不均衡时，本发明方法对核酸序列数据库中每类标签核酸序列使用各自分割步长将每类标签核酸序列分割成等长度的短序列组成短核酸序列均衡样本，避免因数据库核酸序列不平衡导致数据库不可用或需要付出额外努力增加核酸样本的现象，增加了机器学习对核酸序列进行分类的适用范围。

由于某一长度的Kmer的全部数目为4

附图说明

为了更清楚地说明本申请实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1是本发明方法的流程图；

图2是本发明实施例的核酸序列数据库中各亚型所对应的数目示意图，共184个亚型，根据亚型的数目从大到小排序，截取前26种亚型名称及其数目；

图3是通过本发明实施例核酸数据库中各亚型的数目求得每种亚型相较于某一特定值的比率结果图，共184种亚型的比率值，根据比率值从大到小排序，截取前20种亚型的名称及比率值；

图4是将本发明实施例核酸数据库中184种亚型的基因组全长DNA序列分割后，根据短序列的亚型数目排序结果示意图，截取前21种亚型名称及所对应的短序列数目；

图5是本发明实施例使用jellyfish软件提取本发明核酸数据库中的所有20mer形成的txt文件的示意图，截取前13种20mer和数目，‘>49’表示下面20mer(‘AAAAAAAAAAAAAAAAAAAA’)序列在数据库中的数目；

图6是本发明将20mer集合按字母顺序从小到大排列形成20mer字典示意图，不存在的20mer序列标记为，从本本发明核酸数据库中提取出的字典数目共计6942463个不同的20mer序列，截取前25种20mer；

图7显示本发明实施例前20个1000bp DNA序列的数字信号示意图，转化后每个序列的信号长度为981，9284表示为20mer字典中第9824位20mer表示的序列；

图8是本发明实施例模型构建流程图；

图9是本发明实施例的模型训练过程示意图；

图10是本发明实施例的待检测核酸序列样本fasta文件部分展示示意图，‘>1’表示序列的名称，下面为其核酸序列，此图展示了两条待检测序列的名称(‘>1’和‘>14’)和核酸序列；

图11是本发明实施例的待检测核酸序列输入到模型后的分类结果示意图，ID为fasta文件中每条序列的名称，Top-1至Top-5为这条ID序列前5种可能的亚型，括号中的值为此种亚型的可靠性分值，分值越高则表明这条ID序列为此亚型的可靠性越高。

具体实施方式

为了使本领域技术人员更好地理解本申请中的技术方案，下面将结合实施例对本发明作进一步说明，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都应当属于本申请保护的范围。下面结合附图及实施例对本发明作进一步描述。

如图1所示，本发明的目的是利用机器学习对核酸序列进行分类，本实施例以HIV-1型病毒为目标，具体方法如下：

S0：收集带有准确分类标签的核酸序列构建核酸序列数据库：收集带有准确分类标签的HIV基因组DNA数据构建基础数据库，在本实施例中分类标签是HIV各种亚型标签。

本数据库来源于HIV sequence database(https://www.hiv.lanl.gov/components/sequence/HIV/search/search.html)中的12825条完整的HIV-1型DNA序列(截止到2020年9月6日)，然后将其中每种亚型含有的DNA序列数目多于2的基因组DNA序列构建核酸序列数据库，共计184种亚型。

S1：将核酸序列数据库中所有核酸序列分割成等长度的短序列，从而形成短核酸序列均衡样本，该短序列的长度为4bp～核酸序列数据库中最短序列的长度，短核酸序列均衡样本适合于机器学习，短核酸序列均衡样本是指最多分类标签数目的短核酸序列数量N

S1-1：判断核酸序列数据库中所有核酸序列是否均衡。在本发明实施例中，由于带有HIV-1_B亚型标签的序列有5703条，带有HIV-1_AHJU亚型标签的序列有2条(未展示)，所以数据是不均衡的，如图2所示。

S1-2-2：对核酸序列数据库中每类标签核酸序列使用各自分割步长将每类标签核酸序列分割成等长度的短序列，每类标签核酸序列分割步长为分割各类标签核酸序列时后移长度，该后移长度不大于所述短序列的长度，分割核酸序列数据库中所有核酸序列获得短序列样本。

(1)获得核酸序列数据库中带有准确分类标签的核酸序列中每类核酸序列数目与特定值的比率，特定值的选择范围在核酸序列数据库中核酸序列中最多标签数目的核酸序列数量N

(2)根据每类核酸的比率值求得此类核酸的分割步长，分割步长＝短序列长度÷比率值。

在本实施例中短序列长度为1000bp，不同步长的计算方法如下：

(1)通过核酸序列数据库中每种亚型的数目求得每种亚型与特定值的比率。特定值的选择范围在数据库中含有核酸序列最多的标签数目(N

(2)根据每类核酸的比率值求得此类核酸的分割步长，分割步长＝短序列长度÷比率值。例如，HIV-1_B亚型的步长为step＝1000bp÷1.0＝1000bp,HIV-1_01BC亚型的步长为step＝1000bp÷9.94＝100.603≈100bp。

根据每种亚型的步长将核酸序列数据库中对应亚型的全长基因组DNA序列分割成1000bp短DNA序列。分割核酸序列数据库后可得到HIV-1_C亚型含1000bp短DNA序列最多，数目为5066条；HIV-1_AHJU亚型含1000bp短DNA序列最少，数目为2954条(在图中未显示)，如图4所示。

S1-3：判断S1-2-1或者S1-2-2获得短核酸序列的样本是否为均衡并适合于机器学习的样本，如果是，则停止S1-2-1或者S1-2-2，如果不是，则重复S1-2-1或者S1-2-2，直到获得均衡并适合于机器学习的短核酸序列的样本。在本实施例中，由于N

S2-1：使用软件提取核酸序列数据库中所有的Kmer组成集合。在本实施例中，使用jellyfish软件提取数据库种的所有20mer。Kmer的取值范围为：4≤Kmer长度≤短序列长度，在本实施例中短序列长度为1000bp，Kmer选择20mer，如图5所示。

S2-2：构建Kmer字典，将S2-1中的所有Kmer按字母顺序从小到大排序，并且排序好后的Kmer中的第一位置插入表示字典中不存在的Kmer。

在本实施例中，构建20mer字典。将S2-1中的所有20mer形成集合并按字母顺序从小到大排序，考虑到在待检测序列中的20mer可能不存在S2-1中的20mer集合中，因此在排序好后的20mer中的第一位置插入表示字典中不存在的20mer。字典的查询方式为20mer:20mer所在字典中的位置。如图6所示。

S3：使用S2中的Kmer字典将S1中的短核酸序列均衡样本转化成数字信号。在本实施例中，使用20mer字典将S1-3中的均衡样本转化成数字信号，根据样本序列的Kmer在Kmer字典中的位置将样本转化成数字信号，如图7所示。

S4：使用S3中的数字信号训练一个能对带有分类标签的核酸序列进行多分类的神经网络分类器，该神经网络分类器用于识别和判断待检测核酸序列在核酸数据库中核酸分类。

S4-1：搭建一个含有输入层、Embedding层、最大池化层、平均化层、中间合并层、全连接层一、全连接层二、输出层的模型。在S4-1中，输入层维度为1×981，维度计算公式为：短DNA序列长度-Kmer+1，在此实施例中得到的信号长度为1000-20+1＝981；Embedding层维度为6942463×100，其中6942463为构建字典中Kmer的种类；最大池化层为获取embs(981×100)中每一列的最大值，得到一维数组max_pool(1×100)；平均池化层为获取embs中每一列的平均值,得到一维数组avg_pool(1×100)；中间合并层为将max_pool与avg_pool两者通过numpy模块中的concatenate函数进行拼接，得到一维数组input(1×200)；全连接层一含有150个神经元；全连接层二含有150个神经元；输出层含有184个神经元，最后采用softmax激活函数把输出层每个神经元的值转化成概率值。如图8所示。

S4-2：取S3数据集中80％的数据作为训练集，剩余20％的数据作为测试集。在选取时，需要保证训练集和测试集中各亚型的比例近似相同。将训练集送入S4-1中搭建的神经网络中，对网络进行反复训练，并利用测试集不断测试。直到该网络的分类表现达到要求，如图9所示。

S4-3：将上一步训练好的模型分类器保存到磁盘中。对后续待检测核酸序列进行分析。

S5：使用S2中的Kmer字典将待检测序列转化成数字信号。

S5-1：获得待检测HIV核酸序列65条，长度1000bp左右，如图10所示。

S5-2：使用S2-2中的20mer字典将S5-1中的65条HIV核酸序列转化成数字信号。

S6：将S5中的数字信号输入到S4中训练好的神经网络分类器中，对待检测核酸序列进行分类。

S7:展示待检测核酸序列的分类结果。展示结果如图11所示，ID为fasta文件中每条序列的ID，Top-1至Top-5为这条ID序列前5种可能的亚型。括号中的值为此种亚型的可靠性分值，分值越高则表明这条ID序列为此亚型的可靠性越高。

本实施例使用本地计算机进行分析，本地计算机配置为:处理器Intel(R)Core(TM)i5-9400 CPU@2.90GHz,内存8.00GB。神经网络训练好后，9秒内可对6万条核酸序列进行分类，分类结果kappa系数＝0.982。并且由于在输出阶段可以看出检测核酸序列前5种可能的亚型，因此可以根据其概率值的分配评估其是否由多种亚型重组而成，若前两种亚型的概率值均在0.4以上则可判定此序列是由前两种亚型重组而成。

应该理解到披露的本发明不仅仅限于描述的特定的方法、方案和物质，因为这些均可变化。还应理解这里所用的术语仅仅是为了描述特定的实施方式方案的目的，而不是意欲限制本发明的范围，本发明的范围仅受限于所附的权利要求。

本领域的技术人员还将认识到，或者能够确认使用不超过常规实验，在本文中所述的本发明的具体的实施方案的许多等价物。这些等价物也包含在所附的权利要求中。