一种宏基因组重叠群无监督聚类方法及系统技术方案

本发明专利技术公开了一种宏基因组重叠群无监督聚类方法，所述方法是在聚类之前首先将各个样本的reads放在一起组建基因文库，然后通过组装工具将reads组装成contigs，根据四联寡核苷酸的频率以及co

【技术实现步骤摘要】
一种宏基因组重叠群无监督聚类方法及系统


[0001]本专利技术涉及生物信息分析
，尤其涉及一种宏基因组重叠群的无监督聚类方法及系统。

技术介绍

[0002]在宏基因组技术出现之前，人们对于微生物的相关研究主要是通过人工对单一微生物物种进行纯培养。但是自然环境中，绝大多数微生物很难或者不能在培养基上纯培养。随着第二代测序技术的发展，宏基因组技术应运而生。它可以直接从自然环境中获取样本全部微生物的遗传物质，而不需要像传统方法一样在培养基上纯培养。这为科学家研究微生物的群落结构，微生物之间的相互作用以及微生物与环境，疾病之间的关系提供了新的研究思路。第二代测序得到的短的鸟枪宏基因组reads片段可以被reads组装工具组装成更长的基因片段contigs，由于组装工具的局限性，不能组装出完整的基因，而只是零散的基因片段。采用机器学习的方法可以将零散的基因片段分类，从而得到完整的基因，以供后续的物种注释和功能分析。
[0003]现有的宏基因组分类方法一般分为两种，有监督分类和无监督聚类方法。早期，宏基因组组装工具不能很好的对测序得到的reads片段进行组装，大多数分类方法的数据对象是reads片段,由于reads片段很短只有50bp到200bp，携带信息很少，很难对reads进行有效的分类，随着组装工具的准确率上升，可以达到百分之98。越来越多的方法针对组装后contigs片段进行分类，这里我们讨论对组装后的contigs片段分类的方法。有监督分类方法，利用已知的基因作为参考，根据基因序列的同源性以及序列组成的相似性进行分类，由于需要自我构建参考数据库和索引，对计算机的内存和硬盘存储空间要求很高，同时由于环境中有大量未知的物种，无法与参考数据库里的序列进行匹配，所以会有大量无法分类的contigs。正相反，无监督聚类方法可以通过样本中序列本身的组成特点或它们的丰度信息或者结合这两种特点，进行聚类，可以得到未知的物种完整基因，从而发现新的物种。目前主流的方法主要有CONCOCT、Maxbin2.0、MetaBAT、MyCC、COCACOLA、DAS tool以及Metawrap。CONCOCT结合序列的组成信息和co-abundance，将序列特征向量化，然后用PCA方法进行降维，用高斯混合模型结合EM算法进行分类,在简单复杂度数据集中(样本中物种数量50左右)表现很好，但是在复杂环境数据集中表现不行。Maxbin2.0和MetaBAT结合序列组成特征以及co-abundance,并通过预先训练的概率分布模型来计算每个序列到聚类中心点的概率，然后分别用改进后的EM算法和k-medoid算法进行分类，其中Maxbin2.0在中等复杂度数据集上(样本中物种数量300左右)表现很好，但是在高复杂度数据集上(样本中物种数量700左右)重建的高质量基因数量会降低，同时无法应用到超高复杂度数据集(物种数量1000以上)。MetaBAT是专门为复杂数据集设计的算法，可以在超高复杂度数据集上取得很好的效果，缺点是算法参数太多，针对不同的数据集要调参否则无法达到预期的效果。MyCC在多个样本之一中结合了基因组特征，标志基因和可选的重叠群覆盖范围,在低和中等复杂度数据集上表现得很好，但是在复杂度高的数据集表现大大下降。COCACOLA相似性度量
没有使用欧式距离而是使用距离，同时通过稀疏正则化结合硬聚类和软聚类的优点，此外COCACOLA还结合了定制的知识来提高聚类准确率，和大多数方法一样无法在复杂环境数据集中获得很好的表现。没有任何一个方法可以在所有的环境样本中获得很好的表现，于是出现了DAS tool以及Metawrap，它们之中可以添加聚类方法例如CONCOCT、MyCC、Maxbin2.0、MetaBAT，然后将这些聚类算法得到的聚类结果进行去冗余的整合可以得到更多重建的基因。DAS tool以及Metawrap包含的这些聚类方法是同时应用在同一数据集，有些方法在复杂数据集表现不好(例如CONCOCT、MyCC)，同时大多数方法在复杂数据集上重建基因的能力有待提高，所以这些方法聚类结果的整合虽然能提升重建基因的个数以及质量，但是效果却有限。
[0004]目前现有的无监督聚类方法主要问题在于：1.在复杂环境下重建基因能力有待提高。2.聚类过程中菌株个数的选取和真实情况相差很大，菌株个数是无监督聚类方法的关键参数，很大程度上影响算法的性能。3.很难对样本中同一物种但是不同菌株进行区分，有两方面原因一个是使用组装工具对reads进行组装，由于同一物种但是不同菌株之间的相似度很高，组装过程中容易产生嵌合体，第二个原因是复杂环境下菌株数量较多，很难进行有效区分。
[0005]针对这些问题，本专利技术提出了一种宏基因组重叠群无监督聚类方法及系统。

技术实现思路

[0006]本专利技术的目的是针对现有技术的缺陷，提供了宏基因组重叠群无监督聚类方法及系统。
[0007]为了实现以上目的，本专利技术采用以下技术方案：
[0008]一种宏基因组重叠群无监督聚类方法，包括步骤：
[0009]S1.将多个reads通过片段的重叠组装成contigs；
[0010]S2.采用特征向量化方法计算contigs的四联寡核苷酸频率和丰度，得到contigs的四联核苷酸频率和丰度；
[0011]S3.设置K-means算法中的初始物种个数K1，并随机选择K1个contigs作为初始聚类中心；
[0012]S4.根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，并将每个contigs分配给计算得到的概率中概率最大的聚类中心；
[0013]S5.使用每个簇的样本均值更新聚类中心；
[0014]S6.重复执行步骤S4-S5，直到聚类中心不再发生变化时执行步骤S7；
[0015]S7.判断聚类得到的簇中每一个contigs到聚类中心序列的相似概率是否超过最小概率阈值，若是，则得到与contigs相对应的序列集bins，并执行步骤S8；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1；
[0016]S8.使用宏基因组质量检测工具对得到的bins进行质量检测，根据质量检测结果筛选出满足预设阈值的bins，并将筛选出的bins放到集合Q中；将不满足预设阈值的bins混合，记为集合S2；
[0017]S9.将得到的集合S1和集合S2整合为集合S，给定物种数量K，同时初始化K个中心点
medoid；
[0018]S10.计算每个contigs和K个中心点medoid之间的概率，将contigs分配到概率最大的medoid代表的类中；
[0019]S11.重复执行步骤S9-S10，直到所有的medoids不再变化时执行步骤S12；
[0020]S12.判断聚类得到的簇中每一个contigs到medoid的相似概率是否超过最小概率阈值，若是，则返回步骤S8，直到集合Q中不再添加新的bins，则本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种宏基因组重叠群无监督聚类方法，其特征在于，包括步骤：S1.将多个reads通过片段的重叠组装成contigs；S2.采用特征向量化方法计算contigs的四联寡核苷酸频率和丰度，得到contigs的四联核苷酸频率和丰度；S3.设置K-means算法中的初始物种个数K1，并随机选择K1个contigs作为初始聚类中心；S4.根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，并将每个contigs分配给计算得到的概率中概率最大的聚类中心；S5.使用每个簇的样本均值更新聚类中心；S6.重复执行步骤S4-S5，直到聚类中心不再发生变化时执行步骤S7；S7.判断聚类得到的簇中每一个contigs到聚类中心序列的相似概率是否超过最小概率阈值，若是，则得到与contigs相对应的序列集bins，并执行步骤S8；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1；S8.使用宏基因组质量检测工具对得到的bins进行质量检测，根据质量检测结果筛选出满足预设阈值的bins，并将筛选出的bins放到集合Q中；将不满足预设阈值的bins混合，记为集合S2；S9.将得到的集合S1和集合S2整合为集合S，给定物种数量K，同时初始化K个中心点medoid；S10.计算每个contigs和K个中心点medoid之间的概率，将contigs分配到概率最大的medoid代表的类中；S11.重复执行步骤S9-S10，直到所有的medoids不再变化时执行步骤S12；S12.判断聚类得到的簇中每一个contigs到medoid的相似概率是否超过最小概率阈值，若是，则返回步骤S8，直到集合Q中不再添加新的bins，则聚类完成；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1。2.根据权利要求1所述的一种宏基因组重叠群无监督聚类方法，所述步骤S4中根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，具体为：根据贝叶斯公式计算两个contigs来自同一基因的后验概率，表示为：其中，T和R分别表示两个contigs来自不同或者相同菌株的两种情况；D表示两个contigs的四联寡核苷酸频率之间的欧氏距离；根据泊松分布评估序列S和聚类中心点G在宏基因组数据K中相似的概率，表示为：P
cov
(S∈G|cov(G
k
))＝Poisson(cov(S
k
)|cov(G
k
))其中，cov(S
k
)和cov(G
k
)表示序列S和聚类中心序列G在宏基因组数据k中的覆盖度，P
cov
(S
k
∈G
K
|cov(G
k
))为泊松概率密度函数，给定均值λ＝cov(G
k
)；假定所有宏基因组样本是被独立测序的，则序列S和聚类心中序列G之间覆盖度相似概率需要考虑所有宏基因组样本，表示为：
其中，M表示宏基因组样本的数量；每个contigs到各个聚类中心的概率，表示为：其中，P(S∈G)表示每个contigs到各个聚类中心的概率。3.根据权利要求2所述的一种宏基因组重叠群无监督聚类方法，其特征在于，所述步骤S9中给定物种数量K，同时初始化K个中心点medoid，具体为：用reading frames工具将样本中的DNA序列翻译为蛋白质序列，在所有已经被翻译或者未被翻译的序列中，识别34种具有分类信息的标志基因的DNA和蛋白质的隐马尔可夫模型profiles；找出标志基因所对应的DNA序列，并结合FragGeneScan软件预测的基因做去冗余的整合；采用HMMER3软件对整合后的基因进行107个单标记分析，对进行单标记分析的基因进行过滤，得到最短的标记基因序列个数，作为聚类中菌株的数量，以及将最短的标记基因序列作为聚类中心点，并进行聚类中心初始化。4.根据权利要求1所述的一种宏基因组重叠群无监督聚类方法，其特征在于，所述步骤S9中将得到的集合S1和集合S2整合为集合S后还包括删除集合S1和集合S2。5.根据权利要求4所述的一种宏基因组重叠群无监督聚类方法，其特征在于，所述步骤S10中还包括对于第i个类中除对应的medoid
i
之外的其他contigs，按顺序计算当其他contigs成为新的medoid时准则函数的值，寻找准则函数最小时对应的点作为新...

【专利技术属性】
技术研发人员：李小波，姜忠俊，
申请(专利权)人：浙江师范大学，
类型：发明
国别省市：