本公开属于生物信息学分析技术领域，尤其涉及一种微生物宏基因组分箱方法及系统。
背景技术：
本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。
随着测序技术的发展，基于培养与分离的微生物研究方法逐渐被宏基因组学研究方法取代。由于宏基因组是对群落中全部微生物进行测序研究，因此要想对群落的单菌序列进行研究，需要将宏基因组按照菌株或者种系进行分类，在宏基因组学将这种分类方法称为分箱(binning)。
发明人发现，现有的分箱算法对序列特征的利用不够全面，无法对微生物序列特征进行深度挖掘，导致分箱精度不高，同时，现有的分箱算法在对序列特征的降维处理中丢失了大量的有用特征，导致降维后的序列特征无法对微生物进行精确描述。
技术实现要素：
本公开为了解决上述问题，提供了一种微生物宏基因组分箱方法及系统，所述方案采用了多种特征融合的方式，对宏基因组的序列特征进行深度挖掘，实现了对宏基因组序列的精确描述，同时，为了提高算法的处理效率，利用vae-gan神经网络对所提取的特征进行特征降维，所述降维方法在降低特征维度的同时，充分保留了序列特征中的必要成分，良好的平衡了分箱精度与分箱时间的关系，并通过dbscan聚类方法，获得更加精确的分箱结果。
根据本公开实施例的第一个方面，提供了一种微生物宏基因组分箱方法，包括：
获取待分箱的微生物宏基因组序列；
对所述宏基因组序列中每条序列进行特征提取，其中，所提取的特征包括四核苷酸频率、rpkm丰度以及kmer覆盖度特征；
将提取的特征输入vae-gan神经网络中进行训练，通过训练将提取的特征编码到vae隐含向量中；
基于所述vae隐含向量中的均值变量，对所述宏基因组序列进行聚类，根据聚类结果实现宏基因组的分箱。
进一步的，所述vae-gan神经网络中的vae网络，具体包括编码层和解码层，所述编码层首先包括一个全连接网络层，并采用relu作为激活函数，该层的输出作为另外两个全连接网络层的输入，两个全连接层的输出分别对应于vae隐含向量中的均值变量和方差变量。
进一步的，所述vae-gan神经网络中的gan判别器网络，具体包括三个网络组成，第一个网络为两层卷积网络；第二个网络为两层全连接网络，其中第一层的激活函数为leakyrelu函数，第二层的激活函数为sigmoid函数；第三个网络为一个卷积网络。
进一步的，所述四核苷酸频率特征的提取，具体包括：对于每条序列，统计其不包含模糊碱基的四核苷酸频率，获得103维的特征向量。
进一步的，所述丰度特征的计算，具体包括计算映射到每条序列的单条测序读长的数量，若果测序读长被映射到n条序列，则每条序列映射读长计数为1/n，通过获得的序列长度和映射读长的总数，对计数结果进行标准化处理，获得rpkm丰度。
根据本公开实施例的第二个方面，提供了一种微生物宏基因组分箱系统，包括：
数据获取单元，其用于获取待分箱的微生物宏基因组序列；
特征提取单元，其用于对所述宏基因组序列中每条序列进行特征提取，其中，所提取的特征包括四核苷酸频率、rpkm丰度以及kmer覆盖度特征；
特征降维单元，其用于将提取的特征输入vae-gan神经网络中进行训练，通过训练将提取的特征编码到vae隐含向量中；
聚类单元，其用于基于所述vae隐含向量中的均值变量，对所述宏基因组序列进行聚类，根据聚类结果实现宏基因组的分箱。
根据本公开实施例的第三个方面，提供了一种电子设备，包括存储器、处理器及存储在存储器上运行的计算机程序，所述处理器执行所述程序时实现所述的一种微生物宏基因组分箱方法。
根据本公开实施例的第四个方面，提供了一种非暂态计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现所述的一种微生物宏基因组分箱方法。
与现有技术相比，本公开的有益效果是：
(1)本公开所述方案提出了一种针对宏基因组分箱方法，所述方案采用了多种特征融合的方式，对宏基因组的序列特征进行深度挖掘，实现了对宏基因组序列的精确描述，同时，为了提高算法的处理效率，利用vae-gan神经网络对所提取的特征进行特征降维，所述降维方法在降低特征维度的同时，充分保留了序列特征中的必要成分，良好的平衡了分箱精度与分箱时间的关系。
(2)本公开所述方案使用dbscan聚类方法，拥有更好的聚类效果，使得分箱结果更精确。
本公开附加方面的优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本公开的实践了解到。
附图说明
构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。
图1为本公开实施例一中所述的一种微生物宏基因组分箱方法流程图。
具体实施方式
下面结合附图与实施例对本公开做进一步说明。
应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
在不冲突的情况下，本公开中的实施例及实施例中的特征可以相互组合。
实施例一：
本实施例的目的是提供一种微生物宏基因组分箱方法。
一种微生物宏基因组分箱方法，包括：
获取待分箱的微生物宏基因组序列；
对所述宏基因组序列中每条序列进行特征提取，其中，所提取的特征包括四核苷酸频率、rpkm丰度以及kmer覆盖度特征；
将提取的特征输入vae-gan神经网络中进行训练，通过训练将提取的特征编码到vae隐含向量中；
基于所述vae隐含向量中的均值变量，对所述宏基因组序列进行聚类，根据聚类结果实现宏基因组的分箱。
具体的，为了便于理解，以下结合附图对本公开所述方案进行详细说明：
如图1所示展示了一种微生物宏基因组分箱方法流程图，所述方案首先计算样本的tnf频率、丰度以及kmer覆盖度等特征，然后将特征输入到vae-gan神经网络中进行训练，经过训练后，上述特征被编码到该神经网络的vae隐含向量，然后使用dbscan算法对隐含向量进行聚类，并使用向量的余弦距离作为聚类距离。之后得到的聚类结果变为分箱结果。
进一步的，所述特征提取过程具体如下：
对于每条序列，统计不包含模糊碱基的四核苷酸频率(tnf)，最终每条序列可以得到一个103维的tnf向量。对于n个序列，结果则为n×103的张量。为了确定丰度，计算映射到每条序列的单条测序读长的数量，如果测序读长被映射到n条序列，则每条序列映射读长计数为1/n。使用序列长度和映射读长的总数对上述计数进行标准化，从而得到rpkm(readsperkilobasesequencepermillionmappedreads)丰度。若为s个样本每个样品包含n条序列，则丰度输出是表的维度为n×s，对每一行的丰度值以该行的丰度之和对丰度进行归一化，模拟通过将softmax应用于丰度输出神经元。
使用dsk工具来计算整个数据集的kmer覆盖度，并将kmer及其相应的计数存储为以(ki,coverage(ki))的形式，coverage(ki)是整个数据集中的kmerki的计数。对于每条序列，检索其中所有kmer的覆盖度。最后，通过z-scaling每个四核苷酸序列来标准化tnf与coverage(kmer)，以增加相对序列间的方差。
进一步的，所述特征降维的过程具体包括：
将上述特征计算结果作为vae-gan网络的输入(其中，所述四核苷酸频率、rpkm丰度以及kmer覆盖度特征通过特征拼接，获得组合的特征向量)，网络首先通过vae网络进行编码，返回得到参数化正态分布生成的数据与基于隐含变量(均值与方差对数)计算的kl散度，之后将上述参数化正态分布生成的数据作为vae解码器的输入进行解码，生成一个噪声向量，使用噪声向量用vae解码器生成噪声数据。接下来gan判别器网络对真实样本特征、使用真实样本特征解码生成的数据以及噪声向量解码生成的噪声数据进行评分并记录评分结果。然后根据以上三个评分计算损失，并对gan判别器网络的参数进行优化。接下来是对vae网络的参数优化，使用vae网络对真实样本特征进行编码解码样本生成，记录返回的生成样本以及隐含变量的kl散度。生成噪声向量，使用vae对噪声向量解码生成噪声数据，使用gan判别器网络对真实样本特征、使用真实样本特征解码生成的数据以及噪声向量解码生成的噪声数据进行评分并记录评分结果，之后综合使用该评分结果隐含变量kl散度计算vae网络损失，并对网络的参数进行优化。上述过程反复迭代直至序列特征处理完成为止。
其中，vae网络的编码层的实现为首先一个全连接网络，并且使用relu函数作为激活函数，该层的输出作为另外两个全连接网络层的输入，这两个全连接层的输出即为隐含的均值与方差对数。之后利用这两个参数进行正态分布生成，这将作为解码层的输入，同时该网络还计算均值与方差对数的kl散度。解码层为包括两个网络，第一个为两层全连接网络，激活函数为relu函数，并使用批量归一化函数。第二个为两层反卷积网络，第一层的激活函数为relu函数，第二层的激活函数为tanh函数。relu与tanh函数如下所示，激活函数保证了神经网络对非线性结构有好的拟合能力。
进一步的，所述relu函数具体如下：
函数解析式：relu(x)＝max(0，x)；
导数函数解析式：
进一步的，所述tanh函数具体如下：
函数解析式：
导数函数解析式：tanh′(x)＝(1-tanh(x))(1 tanh(x))；
进一步的，所述gan判别器网络由三个网络组成，第一个网络为两层卷积网络，使用leakyrelu函数并设置阈值为0.01。第二个网络为两层全连接网络，其中第一层的激活函数为leakyrelu函数，第二层的激活函数为sigmoid函数。第三个网络为一个卷积网络。与relu函数不同的是leakyrelu函数只有当输入超过阈值时，该函数才会被激活。sigmoid函数如下所示：
函数解析式：
导数函数解析式：f′(x)＝f(x)(1-f(x))
进一步的，所述损失函数具体如下：
gan判别器网络的损失函数如下：
loss_gan＝0.5×(score-real-1)2.mean() 0.5×(score-fake0)2.mean() 0.5×(score_fake1)2.mean()
另外定义损失函数loss_g与loss_gd如下：
loss_g＝0.5×(0.01×(data_fake-data_real)2.sum() (score_fake-score_real)2)/batchsize
loss_gd＝0.5×(score_fake-1)2.mean()
vae网络的损失函数为：0.01×loss_g loss_gd。
进一步的，所述vae-gan神经网络的优化算法使用adam方法进行参数优化。
进一步的，所述聚类过程具体包括：在完成上述神经网络的训练，获得特征降维结果后，使用dbscan聚类方法对网络中的隐含变量(均值与方差对数)中的均值进行聚类，从而实现序列的分箱，其中，所述聚类方法采用的距离度量使用余弦距离。
实施例二：
本实施例的目的是一种微生物宏基因组分箱系统，包括：
数据获取单元，其用于获取待分箱的微生物宏基因组序列；
特征提取单元，其用于对所述宏基因组序列中每条序列进行特征提取，其中，所提取的特征包括四核苷酸频率、rpkm丰度以及kmer覆盖度特征；
特征降维单元，其用于将提取的特征输入vae-gan神经网络中进行训练，通过训练将提取的特征编码到vae隐含向量中；
聚类单元，其用于基于所述vae隐含向量中的均值变量，对所述宏基因组序列进行聚类，根据聚类结果实现宏基因组的分箱。
在更多实施例中，还提供：
一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成实施例一中所述的方法。为了简洁，在此不再赘述。
应理解，本实施例中，处理器可以是中央处理单元cpu，处理器还可以是其他通用处理器、数字信号处理器dsp、专用集成电路asic，现成可编程门阵列fpga或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
存储器可以包括只读存储器和随机存取存储器，并向处理器提供指令和数据、存储器的一部分还可以包括非易失性随机存储器。例如，存储器还可以存储设备类型的信息。
一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成实施例一中所述的方法。
实施例一中的方法可以直接体现为硬件处理器执行完成，或者用处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器、闪存、只读存储器、可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器，处理器读取存储器中的信息，结合其硬件完成上述方法的步骤。为避免重复，这里不再详细描述。
本领域普通技术人员可以意识到，结合本实施例描述的各示例的单元即算法步骤，能够以电子硬件或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
上述实施例提供的一种微生物宏基因组分箱方法及系统可以实现，具有广阔的应用前景。
以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。
上述虽然结合附图对本公开的具体实施方式进行了描述，但并非对本公开保护范围的限制，所属领域技术人员应该明白，在本公开的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本公开的保护范围以内。
1.一种微生物宏基因组分箱方法，其特征在于，包括：
获取待分箱的微生物宏基因组序列；
对所述宏基因组序列中每条序列进行特征提取，其中，所提取的特征包括四核苷酸频率、rpkm丰度以及kmer覆盖度特征；
将提取的特征输入vae-gan神经网络中进行训练，通过训练将提取的特征编码到vae隐含向量中；
基于所述vae隐含向量中的均值变量，对所述宏基因组序列进行聚类，根据聚类结果实现宏基因组的分箱。
2.如权利要求1所述的一种微生物宏基因组分箱方法，其特征在于，所述vae-gan神经网络中的vae网络，具体包括编码层和解码层，所述编码层首先包括一个全连接网络层，并采用relu作为激活函数，该层的输出作为另外两个全连接网络层的输入，两个全连接层的输出分别对应于vae隐含向量中的均值变量和方差变量。
3.如权利要求1所述的一种微生物宏基因组分箱方法，其特征在于，所述vae-gan神经网络中的gan判别器网络，具体包括三个网络组成，第一个网络为两层卷积网络；第二个网络为两层全连接网络，其中第一层的激活函数为leakyrelu函数，第二层的激活函数为sigmoid函数；第三个网络为一个卷积网络。
4.如权利要求1所述的一种微生物宏基因组分箱方法，其特征在于，所述四核苷酸频率特征的提取，具体包括：对于每条序列，统计其不包含模糊碱基的四核苷酸频率，获得103维的特征向量。
5.如权利要求1所述的一种微生物宏基因组分箱方法，其特征在于，所述丰度特征的计算，具体包括计算映射到每条序列的单条测序读长的数量，若果测序读长被映射到n条序列，则每条序列映射读长计数为1/n，通过获得的序列长度和映射读长的总数，对计数结果进行标准化处理，获得rpkm丰度。
6.如权利要求1所述的一种微生物宏基因组分箱方法，其特征在于，所述基于所述vae隐含向量中的均值变量，对所述宏基因组序列进行聚类，所述聚类方法采用dbscan算法。
7.如权利要求6所述的一种微生物宏基因组分箱方法，其特征在于，所述dbscan算法中采用的距离度量为余弦距离。
8.一种微生物宏基因组分箱系统，其特征在于，包括：
数据获取单元，其用于获取待分箱的微生物宏基因组序列；
特征提取单元，其用于对所述宏基因组序列中每条序列进行特征提取，其中，所提取的特征包括四核苷酸频率、rpkm丰度以及kmer覆盖度特征；
特征降维单元，其用于将提取的特征输入vae-gan神经网络中进行训练，通过训练将提取的特征编码到vae隐含向量中；
聚类单元，其用于基于所述vae隐含向量中的均值变量，对所述宏基因组序列进行聚类，根据聚类结果实现宏基因组的分箱。
9.一种电子设备，包括存储器、处理器及存储在存储器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1-7任一项所述的一种微生物宏基因组分箱方法。
10.一种非暂态计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1-7任一项所述的一种微生物宏基因组分箱方法。
技术总结