CN112466404B - 一种宏基因组重叠群无监督聚类方法及系统 - Google Patents
一种宏基因组重叠群无监督聚类方法及系统 Download PDFInfo
- Publication number
- CN112466404B CN112466404B CN202011474316.6A CN202011474316A CN112466404B CN 112466404 B CN112466404 B CN 112466404B CN 202011474316 A CN202011474316 A CN 202011474316A CN 112466404 B CN112466404 B CN 112466404B
- Authority
- CN
- China
- Prior art keywords
- probability
- clustering
- sequence
- contigs
- cov
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 72
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 99
- 239000003550 marker Substances 0.000 claims abstract description 42
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 29
- 238000001514 detection method Methods 0.000 claims abstract description 22
- 238000004458 analytical method Methods 0.000 claims abstract description 14
- 230000010354 integration Effects 0.000 claims description 24
- 239000012634 fragment Substances 0.000 claims description 18
- 230000006870 function Effects 0.000 claims description 18
- 239000000203 mixture Substances 0.000 claims description 15
- 102000004169 proteins and genes Human genes 0.000 claims description 13
- 108091028043 Nucleic acid sequence Proteins 0.000 claims description 12
- 238000004364 calculation method Methods 0.000 claims description 12
- 108020004414 DNA Proteins 0.000 claims description 8
- 238000001914 filtration Methods 0.000 claims description 6
- 238000012216 screening Methods 0.000 claims description 6
- 238000012217 deletion Methods 0.000 claims description 5
- 230000037430 deletion Effects 0.000 claims description 5
- 108091034117 Oligonucleotide Proteins 0.000 abstract description 11
- 230000001580 bacterial effect Effects 0.000 abstract description 2
- JLCPHMBAVCMARE-UHFFFAOYSA-N [3-[[3-[[3-[[3-[[3-[[3-[[3-[[3-[[3-[[3-[[3-[[5-(2-amino-6-oxo-1H-purin-9-yl)-3-[[3-[[3-[[3-[[3-[[3-[[5-(2-amino-6-oxo-1H-purin-9-yl)-3-[[5-(2-amino-6-oxo-1H-purin-9-yl)-3-hydroxyoxolan-2-yl]methoxy-hydroxyphosphoryl]oxyoxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(5-methyl-2,4-dioxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(6-aminopurin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(6-aminopurin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(6-aminopurin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(6-aminopurin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxyoxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(5-methyl-2,4-dioxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(4-amino-2-oxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(5-methyl-2,4-dioxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(5-methyl-2,4-dioxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(6-aminopurin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(6-aminopurin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(4-amino-2-oxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(4-amino-2-oxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(4-amino-2-oxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(6-aminopurin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl]oxy-5-(4-amino-2-oxopyrimidin-1-yl)oxolan-2-yl]methyl [5-(6-aminopurin-9-yl)-2-(hydroxymethyl)oxolan-3-yl] hydrogen phosphate Polymers Cc1cn(C2CC(OP(O)(=O)OCC3OC(CC3OP(O)(=O)OCC3OC(CC3O)n3cnc4c3nc(N)[nH]c4=O)n3cnc4c3nc(N)[nH]c4=O)C(COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3COP(O)(=O)OC3CC(OC3CO)n3cnc4c(N)ncnc34)n3ccc(N)nc3=O)n3cnc4c(N)ncnc34)n3ccc(N)nc3=O)n3ccc(N)nc3=O)n3ccc(N)nc3=O)n3cnc4c(N)ncnc34)n3cnc4c(N)ncnc34)n3cc(C)c(=O)[nH]c3=O)n3cc(C)c(=O)[nH]c3=O)n3ccc(N)nc3=O)n3cc(C)c(=O)[nH]c3=O)n3cnc4c3nc(N)[nH]c4=O)n3cnc4c(N)ncnc34)n3cnc4c(N)ncnc34)n3cnc4c(N)ncnc34)n3cnc4c(N)ncnc34)O2)c(=O)[nH]c1=O JLCPHMBAVCMARE-UHFFFAOYSA-N 0.000 abstract 1
- 241000894007 species Species 0.000 description 39
- 230000000694 effects Effects 0.000 description 13
- 235000021443 coca cola Nutrition 0.000 description 11
- 244000005700 microbiome Species 0.000 description 8
- 239000002609 medium Substances 0.000 description 7
- 239000002773 nucleotide Substances 0.000 description 5
- 125000003729 nucleotide group Chemical group 0.000 description 5
- 230000008569 process Effects 0.000 description 4
- 238000012163 sequencing technique Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000011084 recovery Methods 0.000 description 3
- 238000010187 selection method Methods 0.000 description 3
- 230000000295 complement effect Effects 0.000 description 2
- 238000007477 logistic regression Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 108700042296 Archaeal Genes Proteins 0.000 description 1
- 108700003860 Bacterial Genes Proteins 0.000 description 1
- 230000006399 behavior Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000007621 cluster analysis Methods 0.000 description 1
- 238000012258 culturing Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000003631 expected effect Effects 0.000 description 1
- 238000010230 functional analysis Methods 0.000 description 1
- 239000001963 growth medium Substances 0.000 description 1
- 238000012165 high-throughput sequencing Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 230000003278 mimic effect Effects 0.000 description 1
- 239000013612 plasmid Substances 0.000 description 1
- 238000013441 quality evaluation Methods 0.000 description 1
- 230000008707 rearrangement Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000003612 virological effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
- G16B40/30—Unsupervised data analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
- G06F18/232—Non-hierarchical techniques
- G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
- G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
- G16B30/20—Sequence assembly
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Artificial Intelligence (AREA)
- General Health & Medical Sciences (AREA)
- Biotechnology (AREA)
- Biophysics (AREA)
- General Engineering & Computer Science (AREA)
- Epidemiology (AREA)
- Databases & Information Systems (AREA)
- Public Health (AREA)
- Software Systems (AREA)
- Bioethics (AREA)
- Probability & Statistics with Applications (AREA)
- General Physics & Mathematics (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开了一种宏基因组重叠群无监督聚类方法，所述方法是在聚类之前首先将各个样本的reads放在一起组建基因文库，然后通过组装工具将reads组装成contigs，根据四联寡核苷酸的频率以及co‑abundance对每个contigs进行特征向量化，然后根据预先训练的概率模型以及递归策略进行聚类，引入CheckM进行聚类结果质量检测，并且用来简化样本的复杂性以及作为算法终止的条件，算法初始化采用marker gene分析得到结果作为菌株数量初始化以及聚类中心序列初始化。
Description
技术领域
本发明涉及生物信息分析技术领域，尤其涉及一种宏基因组重叠群的无监督聚类方法及系统。
背景技术
在宏基因组技术出现之前，人们对于微生物的相关研究主要是通过人工对单一微生物物种进行纯培养。但是自然环境中，绝大多数微生物很难或者不能在培养基上纯培养。随着第二代测序技术的发展，宏基因组技术应运而生。它可以直接从自然环境中获取样本全部微生物的遗传物质，而不需要像传统方法一样在培养基上纯培养。这为科学家研究微生物的群落结构，微生物之间的相互作用以及微生物与环境，疾病之间的关系提供了新的研究思路。第二代测序得到的短的鸟枪宏基因组reads片段可以被reads组装工具组装成更长的基因片段contigs，由于组装工具的局限性，不能组装出完整的基因，而只是零散的基因片段。采用机器学习的方法可以将零散的基因片段分类，从而得到完整的基因，以供后续的物种注释和功能分析。
现有的宏基因组分类方法一般分为两种，有监督分类和无监督聚类方法。早期，宏基因组组装工具不能很好的对测序得到的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)，同时大多数方法在复杂数据集上重建基因的能力有待提高，所以这些方法聚类结果的整合虽然能提升重建基因的个数以及质量，但是效果却有限。
目前现有的无监督聚类方法主要问题在于：1.在复杂环境下重建基因能力有待提高。2.聚类过程中菌株个数的选取和真实情况相差很大，菌株个数是无监督聚类方法的关键参数，很大程度上影响算法的性能。3.很难对样本中同一物种但是不同菌株进行区分，有两方面原因一个是使用组装工具对reads进行组装，由于同一物种但是不同菌株之间的相似度很高，组装过程中容易产生嵌合体，第二个原因是复杂环境下菌株数量较多，很难进行有效区分。
针对这些问题，本发明提出了一种宏基因组重叠群无监督聚类方法及系统。
发明内容
本发明的目的是针对现有技术的缺陷，提供了宏基因组重叠群无监督聚类方法及系统。
为了实现以上目的，本发明采用以下技术方案：
一种宏基因组重叠群无监督聚类方法，包括步骤：
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。
进一步的，所述步骤S4中根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，具体为：
根据贝叶斯公式计算两个contigs来自同一基因的后验概率，表示为：
其中，T和R分别表示两个contigs来自不同或者相同菌株的两种情况；D表示两个contigs的四联寡核苷酸频率之间的欧氏距离；
根据泊松分布评估序列S和聚类中心点G在宏基因组数据K中相似的概率，表示为：
Pcov(S∈G|cov(Gk))＝Poisson(cov(Sk)|cov(Gk))
其中，cov(Sk)和cov(Gk)表示序列S和聚类中心序列G在宏基因组数据k中的覆盖度，Pcov(Sk∈GK|cov(Gk))为泊松概率密度函数，给定均值λ＝cov(Gk)；
假定所有宏基因组样本是被独立测序的，则序列S和聚类心中序列G之间覆盖度相似概率需要考虑所有宏基因组样本，表示为：
其中，M表示宏基因组样本的数量；
每个contigs到各个聚类中心的概率，表示为：
其中，P(S∈G)表示每个contigs到各个聚类中心的概率。
进一步的，所述步骤S9中给定物种数量K，同时初始化K个中心点medoid，具体为：
用reading frames工具将样本中的DNA序列翻译为蛋白质序列，在所有已经被翻译或者未被翻译的序列中，识别34种具有分类信息的标志基因的DNA和蛋白质的隐马尔可夫模型profiles；
找出标志基因所对应的DNA序列，并结合FragGeneScan软件预测的基因做去冗余的整合；
采用HMMER3软件对整合后的基因进行107个单标记分析，对进行单标记分析的基因进行过滤，得到最短的标记基因序列个数，作为聚类中菌株的数量，以及将最短的标记基因序列作为聚类中心点，并进行聚类中心初始化。
进一步的，所述步骤S9中将得到的集合S1和集合S2整合为集合S后还包括删除集合S1和集合S2。
进一步的，所述步骤S10中还包括对于第i个类中除对应的medoidi之外的其他contigs，按顺序计算当其他contigs成为新的medoid时准则函数的值，寻找准则函数最小时对应的点作为新的medoid。
相应的，还提供一种宏基因组重叠群无监督聚类系统，包括：
组装模块，用于将多个reads通过片段的重叠组装成contigs；
第一计算模块，用于采用特征向量化方法计算contigs的四联寡核苷酸频率和丰度，得到contigs的四联核苷酸频率和丰度；
选择模块，用于设置K-means算法中的初始物种个数K1，并随机选择K1个contigs作为初始聚类中心；
第二计算模块，用于根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，并将每个contigs分配给计算得到的概率中概率最大的聚类中心；
更新模块，用于使用每个簇的样本均值更新聚类中心；
第一判断模块，用于判断聚类得到的簇中每一个contigs到聚类中心序列的相似概率是否超过最小概率阈值，若是，则得到与contigs相对应的序列集bins；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1；
检测模块，用于使用宏基因组质量检测工具对得到的bins进行质量检测，根据质量检测结果筛选出满足预设阈值的bins，并将筛选出的bins放到集合Q中；将不满足预设阈值的bins混合，记为集合S2；
整合模块，用于将得到的集合S1和集合S2整合为集合S，给定物种数量K，同时初始化K个中心点medoid；
第三计算模块，用于计算每个contigs和K个中心点medoid之间的概率，将contigs分配到概率最大的medoid代表的类中；
第二判断模块，用于判断聚类得到的簇中每一个contigs到medoid的相似概率是否超过最小概率阈值，若是，则返回检测模块，直到集合Q中不再添加新的bins，则聚类完成；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1。
进一步的，所述第二计算模块中根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，具体为：
根据贝叶斯公式计算两个contigs来自同一基因的后验概率，表示为：
其中，T和R分别表示两个contigs来自不同或者相同菌株的两种情况；D表示两个contigs的四联寡核苷酸频率之间的欧氏距离；
根据泊松分布评估序列S和聚类中心点G在宏基因组数据K中相似的概率，表示为：
Pcov(S∈G|cov(Gk))＝Poisson(cov(Sk)|cov(Gk))
其中，cov(Sk)和cov(Gk)表示序列S和聚类中心序列G在宏基因组数据k中的覆盖度，Pcov(Sk∈GK|cov(Gk))为泊松概率密度函数，给定均值λ＝cov(Gk)；
假定所有宏基因组样本是被独立测序的，则序列S和聚类心中序列G之间覆盖度相似概率需要考虑所有宏基因组样本，表示为：
其中，M表示宏基因组样本的数量；
每个contigs到各个聚类中心的概率，表示为：
其中，P(S∈G)表示每个contig到各个聚类中心的概率。
进一步的，所述整合模块中给定物种数量K，同时初始化K个中心点medoid，具体为：
用reading frames工具将样本中的DNA序列翻译为蛋白质序列，在所有已经被翻译或者未被翻译的序列中，识别34种具有分类信息的标志基因的DNA和蛋白质的隐马尔可夫模型profiles；
找出标志基因所对应的DNA序列，并结合FragGeneScan软件预测的基因做去冗余的整合；
采用HMMER3软件对整合后的基因进行107个单标记分析，对进行单标记分析的基因进行过滤，得到最短的标记基因序列个数，作为聚类中菌株的数量，以及将最短的标记基因序列作为聚类中心点，并进行聚类中心初始化。
进一步的，所述整合模块中将得到的集合S1和集合S2整合为集合S后还包括删除集合S1和集合S2。
进一步的，所述第三计算模块中还包括对于第i个类中除对应的medoidi之外的其他contigs，按顺序计算当其他contigs成为新的medoid时准则函数的值，寻找准则函数最小时对应的点作为新的medoid。
与现有技术相比，本发明提出了一个新的聚类方法，它是基于序列的组成特征四联寡核苷酸的频率以及联合丰度co-abundance进行特征向量化并通过概率模型结合k-means以及k-medoids算法进行聚类，并使用CheckM工具通过递归的策略对样本的复杂性进行简化，从而使聚类达到更好的效果。并提出了新的菌株数量选取的方法。我们的方法相比最先进的五种方法CONCOCT、Maxbin2.0、MetaBAT、MyCC、COCACOLA在CAMI模拟数据集以及在两个真实数据集上都取得很好的效果，能获得更多高质量的基因。尤其是在复杂环境下可以取得更好的效果。
附图说明
图1是实施例一提供的一种宏基因组重叠群无监督聚类方法流程图；
图2是实施例一提供的样本物种数量选取方法流程图；
图3是实施例二提供的本发明于现有技术在CAMI模拟数据集上的结果示意图。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。
本发明的目的是针对现有技术的缺陷，提供了一种宏基因组重叠群无监督聚类方法及系统。
实施例一
本实施例提供一种宏基因组重叠群无监督聚类方法，包括步骤：
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.使用CheckM工具对得到的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。
本实施例在聚类之前需要将各个样本的reads放在一起组建基因文库，然后通过组装工具组装成contigs，根据四联寡核苷酸的频率以及co-abundance对每个contigs进行特征向量化，然后根据预先训练的概率模型以及递归策略进行聚类，引入CheckM进行聚类结果质量检测，并且用来简化样本的复杂性以及作为算法终止的条件。算法初始化我们采用marker gene分析得到物种数量以及初始化序列，整个方法的流程如图1所示。
在步骤S1中，将多个reads通过片段的重叠组装成contigs。
reads：高通量测序平台产生的较短的原始序列称为reads。
contigs：通过拼接软件基于reads之间的overlap区间，拼接获得的序列称为contigs(重叠群)。
当组装成contigs后，采用一种用于聚类的递归策略，它可以不断减少样本的复杂度从而提升从复杂环境下重建基因的能力，一种用于聚类的递归策略主要分为以下两个阶段，第一阶段采用K-means算法并结合提出的新的K值选取方法来进行聚类。然后加入质量评价工具CheckM筛选出满足阈值的bins，以此来减低环境样本的复杂度，再将剩下的bins混合进行第二阶段的聚类；第一阶段具体如步骤S2-S7。第二阶段采用K-medoids算法进行聚类，对第一阶段剩下的contigs通过新方法的K值选取和初始聚类中心序列的选取进行更准确的聚类然后再通过CheckM筛选出满足阈值的bins，然后将剩下的bins混合再进行聚类，如此递归直到CheckM不能筛选出满足阈值的bins；第二阶段具体如步骤S8-S12。
在步骤S2中，采用特征向量化方法计算contigs的四联寡核苷酸频率和丰度，得到contigs的四联核苷酸频率和联合丰度co-abundance。
四联寡核苷酸频率是在任何contigs的正向和反向互补链上使用一个bp滑动窗口扫描所有可能的四聚体(即四个连续的核苷酸)的数目来计算的。具有非核苷酸符号的四聚体(即不是A、T、C或G)被丢弃。由于可以从基因组的任何一条链中获得DNA片段，因此将一个四聚体及其反向互补序列的频率结合起来，总共产生136种可能的四聚体。我们用Zi＝(Zi,1,…,Zi,V)来表示每个contig的组成向量，其中V表示136维度，i的取值范围为1到N，其中N为聚类样本中菌株的数量。联合丰度co-abundance向量仅是映射到contig的M个样本中每个样本中reads的平均碱基数。我们用Yi＝(Yi,1,…,Yi,M)表示每个contig的联合丰度向量，其中M为样本的数量，Yi,n即为第n个样本中的reads映射到contig的平均碱基数量，
在步骤S3中，设置K-means算法中的初始物种个数K1，并随机选择K1个contigs作为初始聚类中心。其中，K1的选择由后面物种数量选取部分得出。
k-means是一种迭代求解的聚类分析算法，预将数据分为K组，则随机选取K个对象作为初始的聚类中心，然后计算每个对象与各个种子聚类中心之间的距离，把每个对象分配给距离它最近的聚类中心。聚类中心以及分配给它们的对象就代表一个聚类。每分配一个样本，聚类的聚类中心会根据聚类中现有的对象被重新计算。这个过程将不断重复直到满足某个终止条件。终止条件可以是没有(或最小数目)对象被重新分配给不同的聚类，没有(或最小数目)聚类中心再发生变化，误差平方和局部最小。
在步骤S4中，根据得到的contigs的四联核苷酸频率和联合丰度，计算K1个contigs中的每个contig到各个聚类中心的概率，并将每个contig分配给计算得到的概率中概率最大的聚类中心。
根据得到的contigs的四联核苷酸频率和联合丰度，计算K1个contigs中的每个contig到各个聚类中心的概率，具体为：
四联核苷酸频率被定义为在给定基因序列中连续的四个核苷酸的频率，已经证明对于基因特征和聚类方面具有特定物种模式。从IMG网站下载3181个细菌和古细菌基因，模拟基因序列片段contigs，然后分别计算从intra-genome(序列来自于同一物种)以及从inter-genome(序列取样来自不同基因)序列提取的四联寡核苷酸频率之间的欧氏距离100万次，从而得到来自同一基因和来自不同基因的四联寡核苷酸频率之间欧氏距离的先验概率分布。根据贝叶斯公式计算两个contigs来自同一基因的后验概率，表示为：
其中，T和R分别表示两个contigs来自不同或者相同菌株的两种情况；D表示两个contigs的四联寡核苷酸频率之间的欧氏距离；设置P(T)＝10*P(R)。
同时不同长短的contigs之间的后验概率可以利用逻辑回归进行近似。
其中Dij代表contigi和contigj四联寡核苷酸频率之间的欧氏距离。b和c是两个逻辑回归的参数，它们由实验数据来估计。
鸟枪测序已经被证明遵循Lander-Waterman模型，它可以利用泊松分布来计算核苷酸的覆盖度，因此根据泊松分布评估序列S和聚类中心点G在宏基因组数据K中相似的概率，表示为：
Pcov(S∈G|cov(Gk))＝Poisson(cov(Sk)|cov(Gk))
其中，cov(Sk)和cov(Gk)表示序列S和聚类中心序列G在宏基因组数据k中的覆盖度，Pcov(Sk∈GK|cov(Gk))为泊松概率密度函数，给定均值λ＝cov(Gk)；
假定所有宏基因组样本是被独立测序的，则序列S和聚类心中序列G之间覆盖度相似概率需要考虑所有宏基因组样本，表示为：
其中，M表示宏基因组样本的数量；
结合Pdist以及pcov为衡量每个contigs到各个聚类中心的概率，表示为：
其中，P(S∈G)表示每个contigs到各个聚类中心的概率。
P(S∈G)公式即为概率模型，概率模型可以根据组装好的contigs片段的四联核苷酸频率以及co-abundance，计算每个contigs序列到各个聚类中心的概率。
在步骤S9中，将得到的集合S1和集合S2整合为集合S，给定物种数量K，同时初始化K个中心点medoid。
整合集合S1和集合S2为集合S，同时清空集合S1和集合S2，并且给定物种数量K(K的选取由后面物种数量选取部分得出)同时初始化K个medoid
样本物种数量选取方法具体为：用六种reading frames工具把样本中的DNA序列翻译为蛋白质序列，在所有已经被翻译或者未被翻译的序列中，识别34种具有分类信息的标志基因的DNA和蛋白质的隐马尔可夫模型profiles。找出对应的标志基因所对应的DNA序列并结合FragGeneScan软件预测的基因，做一个去冗余的整合，然后用HMMER3软件对整合后的基因进行107个single-copy marker分析，然后对这些基因进行过滤，得到最短的标记基因序列个数作为聚类中菌株的数量，以及用这些序列作为聚类中心点，如图2所示。从而得到步骤S3中的K1参数和步骤S9中的K参数，并对第二阶段的算法进行聚类中心初始化，最终物种个数由算法递归得出。
在步骤S10中，计算每个contig和K个中心点medoid之间的概率，将contigs分配到概率最大的medoid代表的类中。
对于第i个类中除对应的medoidi之外的所有其他contigs，按顺序计算当其成为新的medoid时准则函数的值，遍历所有可能，寻找准则函数最小时对应的点作为新的medoid。
在步骤S12中，判断聚类得到的簇中每一个contig到medoid的相似概率是否超过最小概率阈值，若是，则返回步骤S8，直到集合Q中不再添加新的bins，则聚类完成；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1。
当聚类完成之后，集合Q中bins的数量既为聚类得到的物种数量，同时也可以得到了每个bin的质量信息。
本实施例相较于现有技术中的CONCOCT、Maxbin2.0以及MetaBAT、MyCC、COCACOLA在模拟数据集以及真实数据集上。结果证明了本实施例提出的方法提高了复杂环境下重建基因的个数和质量。
实施例二
本实施例提供一种宏基因组重叠群无监督聚类方法与实施例一的不同之处在于：
本实施例对本发明的方法在三种不同复杂度的数据集上进行测试，同时与其他最先进的聚类工具CONCOCT、Maxbin2.0、MetaBAT、MyCC、COCACOLA在同一数据集上进行比较。
数据采用CAMI论文里的专门用于基准测试的模拟数据集，它是为了对各个聚类算法有一个统一的评价标准而产生的。其中分为低复杂度数据集(40个基因组和20个环状元件)，中等复杂度数据集(132个基因组和100个环状元件)，以及高复杂度数据集(596个基因组和478个环状元件)，这些数据集为从最新测序的大约700个微生物分离物和600个不同于公共基因组所代表的菌株、物种、属或顺序的环状元件基因组中。同时它们充分模拟了真实环境下的情形，其中包括大量模切相关的菌株，质粒DNA以及病毒序列等等。
在三种不同复杂度的数据集上，同时与其他最先进的聚类工具CONCOCT、Maxbin2.0、MetaBAT、MyCC、COCACOLA在同一数据集上进行比较。
首先统计准确率大于90％、召回率大于30％的bins数量(满足这个条件被认为是质量好的bins，一般认为准确率大于90％，召回率大于30％的bins是同一物种下有不同的菌株的情况)。由于短contigs上的组成特征不明显会影响聚类效果，因此针对1500bp的contigs进行聚类，长度低于1500bp的contigs不参与聚类。同时设置聚类结束后每个簇中所有contigs到聚类中心的最小概率阈值为80％，以及使用CheckM工具时的预设阈值为准确率均为大于90％，召回率为90％，60％，30％。结果如图3所示。在这些聚类方法中，可以看出本发明的方法在几乎每个召回率阈值中得到最多的基因数量，无论是在中等复杂度数据集，或者是高复杂度数据集。在低复杂度数据集中，在召回率90％以上CONCOCT比我们的方法表现得更好，这可能是因为k-means算法在低复杂度数据集中受脏数据影响比较大，所以第一阶段聚类效果不好。第二阶段用CheckM检测时筛选出的基因个数不多而导致的。在中等复杂度和高复杂度数据集中都表现得很好，获得了更多高质量的bins。尤其是在高复杂度数据集中，本发明的方法比其他五种方法效果都要更好。算法的结果已经带有CheckM给出的各个bins的质量信息。从表1可以看出我们的方法在准确率大于90％，召回率大于30％时重建的基因个数是最多的，在准确率大于95％，召回率大于95％时重建的基因个数只在低复杂度数据集上比CONCOCT方法少，在其他数据集上重建基因个数都是最多的。
表1
图3为CONCOCT、Maxbin2.0、MetaBAT、MyCC、COCACOLA以及本方法在CAMI模拟数据集上的表现。其中图3(a)、图3(b)、图3(c)分别为六种方法在简单复杂度数据集、中等复杂度数据集、高复杂度数据集上当Precision>＝90％时，各个recall阈值的物种数量。表1为六种聚类方法在三种模拟数据集上当Precision>＝95％并且recall>＝95％时的重建的物种数量以及当Precision>＝90％并且recall>＝30％时重建的物种数量。
与现有技术相比，本实施例基于序列的组成特征四联寡核苷酸的频率以及联合覆盖度co-abundance进行特征向量化并通过概率模型结合k-means以及k-medoids算法进行聚类，并使用CheckM通过递归的策略对样本的复杂性进行简化，从而使聚类达到更好的效果。并提出了新的菌株数量选取的方法。我们的方法相比最先进的五种方法CONCOCT、Maxbin2.0、MetaBAT、MyCC、COCACOLA在CAMI模拟数据集以及在两个真实数据集上都取得很好的效果，能获得更多高质量的基因。尤其是在复杂环境下可以取得更好的效果。
实施例三
本实施例提供一种宏基因组重叠群无监督聚类系统，包括：
组装模块，用于将多个reads通过片段的重叠组装成contigs；
第一计算模块，用于采用特征向量化方法计算contigs的四联寡核苷酸频率和丰度，得到contigs的四联核苷酸频率和联合丰度；
选择模块，用于设置K-means算法中的初始物种个数K1，并随机选择K1个contigs作为初始聚类中心；
第二计算模块，用于根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contig到各个聚类中心的概率，并将每个contig分配给计算得到的概率中概率最大的聚类中心；
更新模块，用于使用每个簇的样本均值更新聚类中心；
第一判断模块，用于判断聚类得到的簇中每一个contigs到聚类中心序列的相似概率是否超过最小概率阈值，若是，则得到与contigs相对应的序列集bins；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1；
检测模块，用于使用CheckM工具对得到的bins进行质量检测，根据质量检测结果筛选出满足预设阈值的bins，并将筛选出的bins放到集合Q中；将不满足预设阈值的bins混合，记为集合S2；
整合模块，用于将得到的集合S1和集合S2整合为集合S，给定物种数量K，同时初始化K个中心点medoid；
第三计算模块，用于计算每个contigs和K个中心点medoid之间的概率，将contigs分配到概率最大的medoid代表的类中；
第二判断模块，用于判断聚类得到的簇中每一个contig到medoid的相似概率是否超过最小概率阈值，若是，则返回检测模块，直到集合Q中不再添加新的bins，则聚类完成；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1。
进一步的，所述第二计算模块中根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，具体为：
根据贝叶斯公式计算两个contigs来自同一基因的后验概率，表示为：
其中，T和R分别表示两个contigs来自不同或者相同菌株的两种情况；D表示两个contigs的四联寡核苷酸频率之间的欧氏距离；
根据泊松分布评估序列S和聚类中心点G在宏基因组数据K中相似的概率，表示为：
Pcov(S∈G|cov(Gk))＝Poisson(cov(Sk)|cov(Gk))
其中，cov(Sk)和cov(Gk)表示序列S和聚类中心序列G在宏基因组数据k中的覆盖度，Pcov(Sk∈GK|cov(Gk))为泊松概率密度函数，给定均值λ＝cov(Gk)；
假定所有宏基因组样本是被独立测序的，则序列S和聚类心中序列G之间覆盖度相似概率需要考虑所有宏基因组样本，表示为：
其中，M表示宏基因组样本的数量；
每个contigs到各个聚类中心的概率，表示为：
其中，P(S∈G)表示每个contigs到各个聚类中心的概率。
进一步的，所述整合模块中给定物种数量K，同时初始化K个中心点medoid，具体为：
用reading frames工具将样本中的DNA序列翻译为蛋白质序列，在所有已经被翻译或者未被翻译的序列中，识别34种具有分类信息的标志基因的DNA和蛋白质的隐马尔可夫模型profiles；
找出标志基因所对应的DNA序列，并结合FragGeneScan软件预测的基因做去冗余的整合；
采用HMMER3软件对整合后的基因进行107个单标记分析，对进行单标记分析的基因进行过滤，得到最短的标记基因序列个数，作为聚类中菌株的数量，以及将最短的标记基因序列作为聚类中心点，并进行聚类中心初始化。
进一步的，所述整合模块中将得到的集合S1和集合S2整合为集合S后还包括删除集合S1和集合S2。
进一步的，所述第三计算模块中还包括对于第i个类中除对应的medoidi之外的其他contigs，按顺序计算当其他contigs成为新的medoid时准则函数的值，寻找准则函数最小时对应的点作为新的medoid。
需要说明的是，本实施例提供的一种宏基因组重叠群无监督聚类系统与实施例一类似，在此不多做赘述。
与现有技术相比，本发明提出了一个新的聚类方法，它是基于序列的组成特征四联寡核苷酸的频率以及联合丰度co-abundance进行特征向量化并通过概率模型结合k-means以及k-medoids算法进行聚类，并使用CheckM工具通过递归的策略对样本的复杂性进行简化，从而使聚类达到更好的效果。并提出了新的菌株数量选取的方法。我们的方法相比最先进的五种方法CONCOCT、Maxbin2.0、MetaBAT、MyCC、COCACOLA在CAMI模拟数据集以及在两个真实数据集上都取得很好的效果，能获得更多高质量的基因。尤其是在复杂环境下可以取得更好的效果。
注意，上述仅为本发明的较佳实施例及所运用技术原理。本领域技术人员会理解，本发明不限于这里所述的特定实施例，对本领域技术人员来说能够进行各种明显的变化、重新调整和替代而不会脱离本发明的保护范围。因此，虽然通过以上实施例对本发明进行了较为详细的说明，但是本发明不仅仅限于以上实施例，在不脱离本发明构思的情况下，还可以包括更多其他等效实施例，而本发明的范围由所附的权利要求范围决定。
Claims (8)
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；
所述步骤S4中根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，具体为：
根据贝叶斯公式计算两个contigs来自同一基因的后验概率，表示为：
其中，T表示两个contigs来自不同菌株的情况，R表示两个contigs来自相同菌株的情况；D表示两个contigs的四联核苷酸频率之间的欧氏距离；
根据泊松分布评估序列S和聚类中心点G在宏基因组数据K中相似的概率，表示为：
Pcov(S∈G|cov(Gk))＝Poisson(cov(Sk)|cov(Gk))
其中，cov(Sk)和cov(Gk)表示序列S和聚类中心序列G在宏基因组数据k中的覆盖度，Pcov(Sk∈GK|cov(Gk))为泊松概率密度函数，给定均值λ＝cov(Gk)；
假定所有宏基因组样本是被独立测序的，则序列S和聚类心中序列G之间覆盖度相似概率需要考虑所有宏基因组样本，表示为：
其中，M表示宏基因组样本的数量；
每个contigs到各个聚类中心的概率，表示为：
其中，P(S∈G)表示每个contigs到各个聚类中心的概率。
2.根据权利要求1所述的一种宏基因组重叠群无监督聚类方法，其特征在于，所述步骤S9中给定物种数量K，同时初始化K个中心点medoid，具体为：
用reading frames工具将样本中的DNA序列翻译为蛋白质序列，在所有已经被翻译或者未被翻译的序列中，识别34种具有分类信息的标志基因的DNA和蛋白质的隐马尔可夫模型profiles；
找出标志基因所对应的DNA序列，并结合FragGeneScan软件预测的基因做去冗余的整合；
采用HMMER3软件对整合后的基因进行107个单标记分析，对进行单标记分析的基因进行过滤，得到最短的标记基因序列个数，作为聚类中菌株的数量，以及将最短的标记基因序列作为聚类中心点，并进行聚类中心初始化。
3.根据权利要求1所述的一种宏基因组重叠群无监督聚类方法，其特征在于，所述步骤S9中将得到的集合S1和集合S2整合为集合S后还包括删除集合S1和集合S2。
4.根据权利要求3所述的一种宏基因组重叠群无监督聚类方法，其特征在于，所述步骤S10中还包括对于第i个类中除对应的medoidi之外的其他contigs，按顺序计算当其他contigs成为新的medoid时准则函数的值，寻找准则函数最小时对应的点作为新的medoid。
5.一种宏基因组重叠群无监督聚类系统，其特征在于，包括：
组装模块，用于将多个reads通过片段的重叠组装成contigs；
第一计算模块，用于采用特征向量化方法计算contigs的四联核苷酸频率和丰度，得到contigs的四联核苷酸频率和丰度；
选择模块，用于设置K-means算法中的初始物种个数K1，并随机选择K1个contigs作为初始聚类中心；
第二计算模块，用于根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，并将每个contigs分配给计算得到的概率中概率最大的聚类中心；
更新模块，用于使用每个簇的样本均值更新聚类中心；
第一判断模块，用于判断聚类得到的簇中每一个contigs到聚类中心序列的相似概率是否超过最小概率阈值，若是，则得到与contigs相对应的序列集bins；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1；
检测模块，用于使用宏基因组质量检测工具对得到的bins进行质量检测，根据质量检测结果筛选出满足预设阈值的bins，并将筛选出的bins放到集合Q中；将不满足预设阈值的bins混合，记为集合S2；
整合模块，用于将得到的集合S1和集合S2整合为集合S，给定物种数量K，同时初始化K个中心点medoid；
第三计算模块，用于计算每个contigs和K个中心点medoid之间的概率，将contigs分配到概率最大的medoid代表的类中；
第二判断模块，用于判断聚类得到的簇中每一个contigs到medoid的相似概率是否超过最小概率阈值，若是，则返回检测模块，直到集合Q中不再添加新的bins，则聚类完成；若否，则得到无法分类的序列，并将无法分类的序列集合记为集合S1；
所述第二计算模块中根据得到的contigs的四联核苷酸频率和丰度，计算K1个contigs中的每个contigs到各个聚类中心的概率，具体为：
根据贝叶斯公式计算两个contigs来自同一基因的后验概率，表示为：
其中，T表示两个contigs来自不同菌株的情况，R表示两个contigs来自相同菌株的情况；D表示两个contigs的四联核苷酸频率之间的欧氏距离；
根据泊松分布评估序列S和聚类中心点G在宏基因组数据K中相似的概率，表示为：
Pcov(S∈G|cov(Gk))＝Poisson(cov(Sk)|cov(Gk))
其中，cov(Sk)和cov(Gk)表示序列S和聚类中心序列G在宏基因组数据k中的覆盖度，Pcov(Sk∈GK|cov(Gk))为泊松概率密度函数，给定均值λ＝cov(Gk)；
假定所有宏基因组样本是被独立测序的，则序列S和聚类心中序列G之间覆盖度相似概率需要考虑所有宏基因组样本，表示为：
其中，M表示宏基因组样本的数量；
每个contigs到各个聚类中心的概率，表示为：
其中，P(S∈G)表示每个contigs到各个聚类中心的概率。
6.根据权利要求5所述的一种宏基因组重叠群无监督聚类系统，其特征在于，所述整合模块中给定物种数量K，同时初始化K个中心点medoid，具体为：
用reading frames工具将样本中的DNA序列翻译为蛋白质序列，在所有已经被翻译或者未被翻译的序列中，识别34种具有分类信息的标志基因的DNA和蛋白质的隐马尔可夫模型profiles；
找出标志基因所对应的DNA序列，并结合FragGeneScan软件预测的基因做去冗余的整合；
采用HMMER3软件对整合后的基因进行107个单标记分析，对进行单标记分析的基因进行过滤，得到最短的标记基因序列个数，作为聚类中菌株的数量，以及将最短的标记基因序列作为聚类中心点，并进行聚类中心初始化。
7.根据权利要求5所述的一种宏基因组重叠群无监督聚类系统，其特征在于，所述整合模块中将得到的集合S1和集合S2整合为集合S后还包括删除集合S1和集合S2。
8.根据权利要求7所述的一种宏基因组重叠群无监督聚类系统，其特征在于，所述第三计算模块中还包括对于第i个类中除对应的medoidi之外的其他contigs，按顺序计算当其他contigs成为新的medoid时准则函数的值，寻找准则函数最小时对应的点作为新的medoid。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011474316.6A CN112466404B (zh) | 2020-12-14 | 2020-12-14 | 一种宏基因组重叠群无监督聚类方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011474316.6A CN112466404B (zh) | 2020-12-14 | 2020-12-14 | 一种宏基因组重叠群无监督聚类方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112466404A CN112466404A (zh) | 2021-03-09 |
| CN112466404B true CN112466404B (zh) | 2024-02-02 |
Family
ID=74804289
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011474316.6A Active CN112466404B (zh) | 2020-12-14 | 2020-12-14 | 一种宏基因组重叠群无监督聚类方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112466404B (zh) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113205856B (zh) * | 2021-06-22 | 2022-07-12 | 南开大学 | 一种微生物宏基因组分箱方法及系统 |
| CN113393898B (zh) * | 2021-06-29 | 2024-01-05 | 中国科学院深圳先进技术研究院 | 一种基于自监督学习的宏基因组重叠群分类方法 |
| CN113611359B (zh) * | 2021-08-13 | 2022-08-05 | 江苏先声医学诊断有限公司 | 一种提高宏基因组纳米孔测序数据菌种组装效率的方法 |
| CN114065866B (zh) * | 2021-11-22 | 2024-04-30 | 吉林大学 | 一种基于参考物种标签约束的宏基因组序列深度聚类方法 |
| CN114694755B (zh) * | 2022-03-28 | 2023-01-24 | 中山大学 | 基因组组装方法、装置、设备及存储介质 |
| CN115511884B (zh) * | 2022-11-15 | 2023-03-07 | 江苏惠汕新能源集团有限公司 | 基于计算机视觉的冲孔复合模表面质量检测方法 |
| CN117995283B (zh) * | 2024-04-03 | 2024-07-23 | 吉林大学 | 一种单样本宏基因组聚类方法、系统、终端及存储介质 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102332064A (zh) * | 2011-10-07 | 2012-01-25 | 吉林大学 | 基于基因条形码的生物物种识别方法 |
| EP2626802A2 (en) * | 2012-02-10 | 2013-08-14 | Tata Consultancy Services Limited | Assembly of metagenomic sequences |
| WO2015048595A1 (en) * | 2013-09-27 | 2015-04-02 | Jay Shendure | Methods and systems for large scale scaffolding of genome assemblies |
| CN106055928A (zh) * | 2016-05-29 | 2016-10-26 | 吉林大学 | 一种宏基因组重叠群的分类方法 |
| CN106599618A (zh) * | 2016-12-23 | 2017-04-26 | 吉林大学 | 一种宏基因组重叠群的无监督分类方法 |
| CN108133122A (zh) * | 2016-12-01 | 2018-06-08 | 深圳华大基因股份有限公司 | 基因聚类方法和基于该方法的宏基因组组装方法和装置 |
| CN109273053A (zh) * | 2018-09-27 | 2019-01-25 | 华中科技大学鄂州工业技术研究院 | 一种高通量测序的微生物数据处理方法 |
| WO2019232357A1 (en) * | 2018-05-31 | 2019-12-05 | Arizona Board Of Regents On Behalf Of The University Of Arizona | Methods for comparative metagenomic analysis |
-
2020
- 2020-12-14 CN CN202011474316.6A patent/CN112466404B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102332064A (zh) * | 2011-10-07 | 2012-01-25 | 吉林大学 | 基于基因条形码的生物物种识别方法 |
| EP2626802A2 (en) * | 2012-02-10 | 2013-08-14 | Tata Consultancy Services Limited | Assembly of metagenomic sequences |
| WO2015048595A1 (en) * | 2013-09-27 | 2015-04-02 | Jay Shendure | Methods and systems for large scale scaffolding of genome assemblies |
| CN106055928A (zh) * | 2016-05-29 | 2016-10-26 | 吉林大学 | 一种宏基因组重叠群的分类方法 |
| CN108133122A (zh) * | 2016-12-01 | 2018-06-08 | 深圳华大基因股份有限公司 | 基因聚类方法和基于该方法的宏基因组组装方法和装置 |
| CN106599618A (zh) * | 2016-12-23 | 2017-04-26 | 吉林大学 | 一种宏基因组重叠群的无监督分类方法 |
| WO2019232357A1 (en) * | 2018-05-31 | 2019-12-05 | Arizona Board Of Regents On Behalf Of The University Of Arizona | Methods for comparative metagenomic analysis |
| CN109273053A (zh) * | 2018-09-27 | 2019-01-25 | 华中科技大学鄂州工业技术研究院 | 一种高通量测序的微生物数据处理方法 |
Non-Patent Citations (3)
| Title |
|---|
| MetaCRS: unsupervised clustering of contigs with the recursive strategy of reducing metagenomic dataset’s complexity;Zhongjun jiang等;《BMC Bioinformatics》;1-16 * |
| Unsupervised Binning of Metagenomic Assembled Contigs Using Improved Fuzzy C-Means Method;Yun Liu等;《IEEE/ACM Transactions on Computational Biology and Bioinformatics》;第14卷(第6期);1459-1467 * |
| 新一代测序技术下海量宏基因组数据分类研究;江源;《中国优秀硕士学位论文全文数据库 (信息科技辑)》(第1期);I140-600 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112466404A (zh) | 2021-03-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112466404B (zh) | 一种宏基因组重叠群无监督聚类方法及系统 | |
| Mitra et al. | Multi-objective evolutionary biclustering of gene expression data | |
| EP2511843B1 (en) | Method and system for calling variations in a sample polynucleotide sequence with respect to a reference polynucleotide sequence | |
| Linder et al. | Predicting RNA-seq coverage from DNA sequence as a unifying model of gene regulation | |
| Wang et al. | SCMarker: ab initio marker selection for single cell transcriptome profiling | |
| JP2023551517A (ja) | 人工知能ベースのがん診断及びがん種予測方法｛Ｍｅｔｈｏｄ ｆｏｒ ｄｉａｇｎｏｓｉｎｇ ａｎｄ ｐｒｅｄｉｃｔｉｎｇ ｃａｎｃｅｒ ｔｙｐｅ ｂａｓｅｄ ｏｎ ａｒｔｉｆｉｃｉａｌ ｉｎｔｅｌｌｉｇｅｎｃｅ ｂａｓｅｄ ｏｎ ａｒｔｉｆｉｃｉａｌ ｉｎｔｅｌｌｉｇｅｎｃｅ｝ | |
| CN115631789A (zh) | 一种基于泛基因组的群体联合变异检测方法 | |
| Patel et al. | Multi-Classifier Analysis of Leukemia Gene Expression From Curated Microarray Database (CuMiDa) | |
| Arevalillo et al. | Exploring correlations in gene expression microarray data for maximum predictive–minimum redundancy biomarker selection and classification | |
| CN103348350B (zh) | 核酸信息处理装置及其处理方法 | |
| Zintzaras et al. | Forest classification trees and forest support vector machines algorithms: Demonstration using microarray data | |
| CN114207727A (zh) | 用于从变体识别数据确定起源细胞的系统和方法 | |
| Shehzadi et al. | Intelligent predictor using cancer-related biologically information extraction from cancer transcriptomes | |
| CN111755074B (zh) | 一种酿酒酵母菌中dna复制起点的预测方法 | |
| CN108182347B (zh) | 一种大规模跨平台基因表达数据分类方法 | |
| Ranek et al. | DELVE: feature selection for preserving biological trajectories in single-cell data | |
| Cutigi et al. | A proposal of a graph-based computational method for ranking significant set of related genes in cancer | |
| Kawaguchi et al. | Learning single-cell chromatin accessibility profiles using meta-analytic marker genes | |
| CN111785319A (zh) | 基于差异表达数据的药物重定位方法 | |
| Wong et al. | A gene selection method for microarray data based on risk genes | |
| Blazadonakis et al. | The linear neuron as marker selector and clinical predictor in cancer gene analysis | |
| KR102405732B1 (ko) | 세포 클러스터링 방법 및 장치 | |
| Masárová | Struktura repeatomu u vybraných zástupců rodu Boechera (brukvovité) | |
| Basher et al. | Heterogeneity-Preserving Discriminative Feature Selection for Subtype Discovery | |
| Fahrenberger et al. | GTestimate: Improving relative gene expression estimation in scRNA-seq using the Good-Turing estimator |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |