CN114065866A - 一种基于参考物种标签约束的宏基因组序列深度聚类方法 - Google Patents
一种基于参考物种标签约束的宏基因组序列深度聚类方法 Download PDFInfo
- Publication number
- CN114065866A CN114065866A CN202111389111.2A CN202111389111A CN114065866A CN 114065866 A CN114065866 A CN 114065866A CN 202111389111 A CN202111389111 A CN 202111389111A CN 114065866 A CN114065866 A CN 114065866A
- Authority
- CN
- China
- Prior art keywords
- clustering
- error
- species
- different
- training
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 64
- 238000012549 training Methods 0.000 claims abstract description 61
- 230000000813 microbial effect Effects 0.000 claims abstract description 45
- 230000006870 function Effects 0.000 claims abstract description 25
- 239000013598 vector Substances 0.000 claims abstract description 19
- 238000004364 calculation method Methods 0.000 claims description 31
- 238000012360 testing method Methods 0.000 claims description 14
- ORILYTVJVMAKLC-UHFFFAOYSA-N Adamantane Natural products C1C(C2)CC3CC1CC2C3 ORILYTVJVMAKLC-UHFFFAOYSA-N 0.000 claims description 4
- 230000004913 activation Effects 0.000 claims description 3
- 239000006185 dispersion Substances 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims 1
- 238000011160 research Methods 0.000 abstract description 5
- 238000013135 deep learning Methods 0.000 abstract description 4
- 241000894007 species Species 0.000 description 96
- 238000004088 simulation Methods 0.000 description 11
- 244000005700 microbiome Species 0.000 description 9
- 238000004458 analytical method Methods 0.000 description 7
- 235000021443 coca cola Nutrition 0.000 description 5
- 108091028043 Nucleic acid sequence Proteins 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 239000012634 fragment Substances 0.000 description 3
- 238000005286 illumination Methods 0.000 description 3
- 238000009423 ventilation Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 108020004414 DNA Proteins 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012163 sequencing technique Methods 0.000 description 2
- 241000606125 Bacteroides Species 0.000 description 1
- 241001112696 Clostridia Species 0.000 description 1
- 241000305071 Enterobacterales Species 0.000 description 1
- 241000588724 Escherichia coli Species 0.000 description 1
- 238000003766 bioinformatics method Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 244000005709 gut microbiome Species 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000007481 next generation sequencing Methods 0.000 description 1
- 102000004169 proteins and genes Human genes 0.000 description 1
- 108090000623 proteins and genes Proteins 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/004—Artificial life, i.e. computing arrangements simulating life
- G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/045—Combinations of networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
- G06N3/084—Backpropagation, e.g. using gradient descent
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Molecular Biology (AREA)
- Computational Linguistics (AREA)
- Software Systems (AREA)
- Mathematical Physics (AREA)
- Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Computing Systems (AREA)
- General Health & Medical Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明提供了一种基于参考物种标签约束的宏基因组序列深度聚类方法，设计了基于参考物种标签约束的深度学习预训练模型。本发明建立了基于不同群落的已知物种的预训练数据库，构建预训练数据库时将每条4mer特征向量分为同一物种、相同属不同物种和不同属不同物种三种情况，并分别研究了三种情况下的样本间序列的4mer特征间的关系；建立了预训练模型的标签约束误差函数，并且使用群落已知标签的数据库进行预训练，针对不同的微生物群落构建不同预训练模型；在用户使用时，只需要针对不同的群落加载所需群落的预训练模型，重新加载模型仅仅等待几次微调步骤的迭代即可得到聚类结果。最终，所述聚类方法能够展现非常优秀的聚类性能。
Description
技术领域
本发明涉及生物信息学分析领域，尤其涉及一种基于参考物种标签约束的宏基因组序列深度聚类方法。
背景技术
微生物是地球上种类最大、数量最多、分布最广的生物群。人们对于微生物的研究主要是建立在纯培养的基础上，然而后来发现99％以上的微生物是不可培养的。为了研究不能培养的微生物，一个全新的理念——宏基因组学应运而生。宏基因组学利用新一代测序技术在不经过实验室培养的情况下，能够获取环境中绝大部分的遗传物质。与传统的测序方法不同，宏基因组测序得到的原始数据是大量的、长度较短的、来源于多种微生物的DNA片段。研究人员根据DNA片段之间的重叠关系可以将这些零碎的短片段组装成长度较长的DNA序列，生物信息学中称这种组装之后的DNA序列为重叠群(contigs)。将这些重叠群依据其物种归属进行分类是宏基因组数据分析中十分重要的一步。
然而，受宏基因组中不同物种间的丰度比、基因组长度等因素的影响，组装后属于不同物种的重叠群数量往往相差很多，因此，宏基因组重叠群数据是一种典型的不平衡数据集。如何对其进行有效地分类，是目前的一个研究难点。
因此，现有技术还有待改进。
发明内容
鉴于上述现有技术的不足，本发明的目的在于提供一种于参考物种标签约束的宏基因组序列深度聚类方法，旨在解决现有技术在进行宏基因组DNA序列聚类时，因相同属的临近物种相似度较高而导致的聚类不准确的问题。
本发明的技术方案如下：
本发明提供一种基于参考物种标签约束的宏基因组序列深度聚类方法，其中，包括步骤：
第一步，预训练步骤，包括：
1.1构建基于不同微生物群落的已知物种的预训练数据库；
1.2建立初始化模型；
1.3针对不同的微生物群落构建不同的预训练模型；
第二步，微调步骤，包括：
2.1计算待测微生物群落的数据集重叠群序列4mer频率，并归一化得到微调的输入特征频率Y；
2.2加载待测微生物群落的预训练模型以及参数；
2.3构建微调误差函数；
2.4确定聚类个数；
2.5微调模型；
2.6完成聚类，得到聚类结果，并根据聚类指标衡量聚类结果。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述步骤1.1具体包括步骤：
a)下载不同微生物群落的已知物种的全基因组序列预训练数据集；
b)将每个物种的全基因组序列从随机起始位置截取随机长度的序列；
c)计算步骤b)中截取的每条序列的4mer频率特征，并进行归一化，得到不同微生物群落的宏基因组预训练4mer频率归一化特征X。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述步骤1.2具体包括步骤：
a)建立具有对称结构的自编码器；
b)选取
函数作为激活函数，加入Dropout函数来调节模型参数和样本量之间的关系；
c)设置模型参数；
d)构建预训练误差函数，其计算公式：
ERRORpre＝ERES+κELCN
其中，ERES表示重构误差，ELCN表示标签约束误差，κ表示用于平衡重构误差ERES和标签约束误差ELCN的超参数。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述重构误差ERES的计算公式为：
其中，xi表示编码网络的输入，f(xi)表示编码网络的输出，g(f(xi))表示解码网络的输出，N1表示总样本的个数；
所述标签约束误差ELCN的计算公式为：
其中，Es表示衡量相同物种间的特征向量间的欧式距离，El表示衡量相同属不相同物种间的特征向量间的距离，Ed表示衡量不同属不同物种间的特征向量间的距离，n1、n2、n3为三种误差下累加的次数且满足
β、λ是标签约束相的超参数。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述步骤1.3中，构建预训练模型具体包括步骤：
a)使用初始化后的网络模型以及参数；
b)加载测试集样本归一化特征X，并送入初始化后网络模型，计算重构误差和标签约束误差；
c)应用反向传播自适应矩估计方法，对不同变化的参数以自适应的学习率进行更新；
d)保存预训练的模型以及参数，定义为待测微生物群落的预训练模型。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述步骤2.3具体包括步骤：
利用深度k-mean聚类的方法，将k-mean聚类的误差加入到预训练得到的待测微生物群落的预训练模型中，得到微调误差函数：
其中，ERES表示重构误差，ECLU表示聚类误差，N2表示待测微生物群落样本的个数，η表示用于平衡重构误差ERES和聚类误差ECLU的超参数。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述重构误差ERES计算公式为：
其中，yi表示微调网络测试集的输入，g(f(yi))表示微调网络测试集的输出；
所述聚类误差ECLU计算公式为：
其中，x表示Ci类中的样本点，θi表示Ci类的聚类中心。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述步骤2.4中，确定聚类个数具体包括步骤：
a)定义Si变量为类内数据到簇质心的平均距离，代表簇Ci中各样本的分散程度，其计算公式为：
其中，|Ci|表示簇类Ci中样本点个数，x表示Ci类中的样本点，θi表示Ci类的聚类中心；
b)定义Mij为簇Ci与簇Cj的距离，其计算公式为：
Mij＝||θi-θj||2
其中，Mij表示簇Ci与簇Cj的质心θi和θj的距离；
c)应用DBI指数确定聚类个数k，其计算公式为：
其中，k为聚类个数变量，取值范围为[0.002*N2,0.005*N2]；
d)选择最优的聚类个数k0，其计算公式为：
其中，k为聚类个数变量。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述步骤2.5具体包括步骤：
a)将测试样本送入网络，得到第一次隐藏表示；
b)将a)中得到的隐藏表示作为样本的特征使用k-mean聚类，聚类个数为k0，初始化中心的方法使用k-means++；
c)将b)中得到的聚类结果分批次随机取出，计算重构误差ERES和聚类误差ECLU，并根据重构误差ERES和聚类误差ECLU计算得到微调误差ERRORfine；
d)对每个样本的微调误差ERRORfine使用Adam优化器微调模型参数，得到微调后的模型。
所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其中，所述步骤2.6具体包括步骤：
a)联合样本的覆盖率信息和网络中间层用于聚类，得到聚类结果
b)使用准确率、召回率和ARI来衡量聚类结果，其中，准确率的计算公式为：
召回率的计算公式为：
ARI的计算公式为：
其中，TP和FP分别表示真正属于同一物种的contigs聚在同一类和不同类中的数目，FN和TN分别表示属于不同物种的contigs聚集成同一类和不同类的数目。
有益效果：本发明提供了一种基于参考物种标签约束的宏基因组序列深度聚类方法(Label-constrained deep clustering，LCDC)，所述方法与目前的同类方法有明显不同，其设计了基于参考物种标签约束的深度学习预训练模型。由于所有的生态系统都存在微生物群落，而且不同环境下的微生物群落存在巨大的差异，相同环境下也会因为温度、酸碱度、光照、通气性等细微差距而使得微生物种类存在巨大差异。本发明利用微生物群落差异的特点，对不同环境下的微生物群落进行分析，建立了基于不同群落的已知物种的预训练数据库，构建预训练数据库时将每条4mer特征向量分为同一物种、相同属不同物种和不同属不同物种三种情况，并分别研究了三种情况下的样本间序列的4mer特征间的关系；建立了预训练模型的标签约束误差函数，并且LCDC使用群落已知标签的数据库进行预训练，针对不同的微生物群落构建不同预训练模型；在用户使用时，只需要针对不同的群落加载所需群落的预训练模型，重新加载模型仅仅等待几次微调步骤的迭代即可得到聚类结果。经过多个数据集上进行分析，并对比LCDC与其他多种方法的聚类结果，在真实数据集中，LCDC能够展现非常优秀的聚类性能。
附图说明
图1为本发明实施例中一种基于参考物种标签约束的宏基因组序列深度聚类方法的模型结构示意图。
图2为本发明实施例中三个代表行物种的特征运动趋势。
图3为本发明实施例中构建预训练模型流程图。
图4为本发明实施例中特征组合的效果图。
图5为本发明实施例中101物种的Sharon模拟数据中各方法的聚类结果对比图。
图6为本发明实施例中20物种的Strain模拟数据集中各方法的聚类结果对比图。
具体实施方式
本发明提供一种基于参考物种标签约束的宏基因组序列深度聚类方法，为使本发明的目的、技术方案及效果更加清楚、明确，以下对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明实施例提供一种基于参考物种标签约束的宏基因组序列深度聚类方法，包括步骤：
S100、预训练步骤；
S200、微调步骤。
本发明基于物种标签约束的宏基因组深度聚类方法(Label-constrained deepclustering，LCDC)与目前的同类方法有明显不同，设计了基于参考物种标签约束的深度学习预训练模型。由于所有的生态系统都存在微生物群落，而且不同环境下的微生物群落存在巨大的差异，相同环境下也会因为温度、酸碱度、光照、通气性等细微差距而使得微生物种类存在巨大差异。如图1所示，本发明利用微生物群落差异的特点，对不同环境下的微生物群落进行分析，建立了基于不同群落的已知物种的预训练数据库，构建预训练数据库时将每条4mer特征向量分为同一物种、相同属不同物种和不同属不同物种三种情况，并分别研究了三种情况下的样本间序列的4mer特征间的关系；建立了预训练模型的标签约束误差函数，并且LCDC使用群落已知标签的数据库进行预训练，针对不同的微生物群落构建不同预训练模型；在用户使用时，只需要针对不同的群落加载所需群落的预训练模型，重新加载模型仅仅等待几次微调步骤的迭代即可得到聚类结果。
在一些实施方式中，S100具体包括以下步骤：
S101、构建基于不同微生物群落的已知物种的预训练数据库；
S102、建立初始化模型；
S103、针对不同的微生物群落构建不同的预训练模型。
在一些实施方式中，所述步骤S101中，构建基于不同微生物群落的已知物种的预训练数据库具体包括步骤：
a)下载不同微生物群落的已知物种的全基因组序列预训练数据集；
b)将每个物种的全基因组序列从随机起始位置截取随机长度的序列；
c)计算步骤b)中截取的每条序列的4mer频率特征，并进行归一化，得到不同微生物群落的宏基因组预训练4mer频率归一化特征X。
在一些优选的实施方式中，所述截取随机长度的序列，其长度在5000bp-8000bp之间。
在一些优选的实施方式中，可根据实际需要确定每条宏基因组序列截取随机长度的序列数目，例如截取800条。
在一些实施方式中，所述步骤S102中，建立初始化模型具体包括步骤：
a)建立具有对称结构的自编码器；
b)选取
函数作为激活函数，加入Dropout函数来调节模型参数和样本量之间的关系；
c)设置模型参数；
d)构建预训练误差函数，其计算公式为(1)：
ERRORpre＝ERES+κELCN (1)
其中，ERES表示重构误差，ELCN表示标签约束误差，κ表示用于平衡重构误差ERES和标签约束误差ELCN的超参数。
在一些优选的实施方式中，步骤S102的步骤a)中，具有对称结构的自编码器的参数设置为输入为136维，隐藏层100*80*50*80*100，输出为136维。
在一些优选的实施方式中，步骤S102的步骤b)中，所述Dropout函数的参数设置为0.5，其目的是防止过拟合，调节模型参数和样本量之间的关系。
在一些优选的实施方式中，步骤S102的步骤c)中，所述模型参数可设置为每一批次训练个数n＝200，自适应矩估计(Adaptive moment estimation，Adam)优化器的自适应的学习率lr_rate＝0.0001，训练周期epoch＝300。
在一些实施方式中，步骤S102的步骤c)中，所述重构误差ERES的计算公式为(2)：
其中，xi表示编码网络的输入，f(xi)表示编码网络的输出，g(f(xi))表示解码网络的输出，N1表示总样本的个数。
在一些实施方式中，步骤S102的步骤c)中，所述标签约束误差ELCN的计算公式为(3)：
其中，Es表示衡量相同物种间的特征向量间的欧式距离，El表示衡量相同属不相同物种间的特征向量间的距离，Ed表示衡量不同属不同物种间的特征向量间的距离，n1、n2、n3为三种误差下累加的次数且满足
β、λ是标签约束相的超参数。
如图2所示，A、B、C分别表示样本的三个物种，其中A、B是相同属的两个相近物种，C表示与A、B不同属的物种。x1表示A物种的一个样本，相距A物种其他样本较远，箭头表示本方法标签约束项的目的：相同物种内样本特征向量表示要近，不同物种的样本特征向量表示尽量要远一些。而且，需要着重考虑相同属的不同物种相距较不同属的不同物种的特征向量要更相似，对相同属的不同物种之间的误差设定更大的超参数。依据上述分析将样本中任意两条的物种信息分为3种情况：第一种是两条特征序列来自同一个物种，第二种是两条特征序列来自不同物种的同一个属，第三种是两条特征序列来自不同物种不同属。对样本任意两条特征向量分为上述三种情况，计算三种情况下的约束误差，计算公式为：
上式中，xi、xj为重叠群任意两条序列的4mer的特征向量，用来衡量相同物种间的特征向量间的欧式距离；用Es来衡量相同物种间的特征向量间的欧式距离，用El来衡量相同属不相同物种间的特征向量间的距离，用Ed来衡量不同属不同物种间的特征向量间的距离。n1、n2、n3为三种误差下累加的次数且满足
β、λ是标签约束相的超参数用来调节三者误差，可设置β＝2，λ＝10。将三种情况下的误差项相加即得到样本总的标签约束误差ELCN。
在一些实施方式中，如图3所示，所述步骤S103中，构建预训练模型具体包括步骤：
a)使用初始化后的网络模型以及参数；
b)加载测试集样本归一化特征X，并送入初始化后网络模型，计算重构误差和标签约束误差；
c)应用反向传播自适应矩估计方法，对不同变化的参数以自适应的学习率进行更新；
d)保存预训练的模型以及参数，定义为待测微生物群落的预训练模型。
在一些实施方式中，S200具体包括以下步骤：
S201、计算待测微生物群落的数据集重叠群序列4mer频率，并归一化得到微调的输入特征频率Y；
S202、加载待测微生物群落的预训练模型以及参数；
S203、构建微调误差函数；
S204、确定聚类个数；
S205、微调模型；
S206、完成聚类，得到聚类结果，并根据聚类指标衡量聚类结果。
在一些实施方式中，所述步骤S203中，构建微调误差函数具体包括步骤：
利用深度k-mean聚类的方法，将k-mean聚类的误差加入到预训练得到的待测微生物群落预训练模型中，得到微调误差函数(4)：
其中，ERES表示重构误差，ECLU表示聚类误差，N2表示待测微生物群落样本的个数，η表示用于平衡重构误差ERES和聚类误差ECLU的超参数。
在前述预训练模型的基础上完成聚类，在此模型中利用一种深度k-mean聚类的方法，将k-mean聚类的误差加入到预训练得到的模型中，得到微调误差函数。
在一些优选的实施方式中，η＝10。
在一些实施方式中，所述重构误差ERES计算公式为(5)：
其中，yi表示微调网络测试集的输入，g(f(yi))表示微调网络测试集的输出；
所述聚类误差ECLU为计算公式为(6)：
其中，x表示Ci类中的样本点，θi表示Ci类的聚类中心。
在一些实施方式中，所述步骤S204中，确定聚类个数具体包括步骤：
a)定义Si变量为类内数据到簇质心的平均距离，代表簇Ci中各样本的分散程度，其计算公式为(7)：
其中，|Ci|表示簇类Ci中样本点个数，x表示Ci类中的样本点，θi表示Ci类的聚类中心；
b)定义Mij为簇Ci与簇Cj的距离，其计算公式为(8)：
Mij＝||θi-θj||2 (8)
其中，Mij表示簇Ci与簇Cj的质心θi和θj的距离；
c)应用DBI指数确定聚类个数k，其计算公式为(9)：
其中，k为聚类个数变量，取值范围为[0.002*N2,0.005*N2]；Si的计算公式为公式(7)；Sj的计算公式与Si的计算公式一致，仅替换为各样本j相应的参数，且j≠i；Mij由公式(8)计算。
d)选择最优的聚类个数k0，其计算公式为(10)：
使用DBI指数确定聚类个数k，将k个类取平均即得到DBI指数，在上述k的取值区间内选择具有最小DBI指数即为k0。
在一些实施方式中，所述步骤S205中，微调模型具体包括步骤：
a)将测试集样本送入网络，得到第一次隐藏表示；
b)将a)中得到的隐藏表示作为样本的特征使用k-mean聚类，聚类个数为k0，初始化中心的方法使用k-means++；
c)将b)中得到的聚类结果分批次随机取出，计算重构误差ERES和聚类误差ECLU，并根据重构误差ERES和聚类误差ECLU计算得到微调误差ERRORfine
d)对每个样本的微调误差ERRORfine使用Adam优化器微调模型参数，得到微调后的模型。
在一些实施方式中，步骤S205的步骤c)中，重构误差ERES和聚类误差ECLU分别取平均再乘以超参数后相加，即得到所述微调误差ERRORfine。
在一些优选的实施方式中，步骤S205的步骤c)中，聚类结果分批次随机取出时，每批设置200个样本。
在一些实施方式中，所述步骤S206中，完成聚类，并根据聚类指标衡量聚类结果具体包括步骤：
a)联合样本的覆盖率信息和网络中间层用于聚类，得到聚类结果
b)使用准确率、召回率和ARI来衡量聚类结果，其中，准确率的计算公式为(11)：
召回率的计算公式为(12)：
ARI的计算公式为(13)：
其中，TP和FP分别表示真正属于同一物种的contigs聚在同一类和不同类中的数目，FN和TN分别表示属于不同物种的contigs聚集成同一类和不同类的数目。
在一些实施方式中，步骤S206的步骤a)中，联合覆盖率信息如图4所示，其中4mer频率特征是136维的，4mer频率特征经过编码网络得到的隐藏层，例如人体肠道微生物预训练网络模型中隐藏层为(100*80*50)，则u的值为50。Species数据集是由96个样本组成的重叠群数据集，因此v的值为96。样本的覆盖率信息和网络中间层联合用于聚类，得到聚类结果
在一些实施方式中，步骤S206的步骤b)中，使用准确率、召回率和调整兰德尔指数(Adjusted Rand index，ARI)来衡量聚类结果。评估标准的重点是同一物种的contigs要聚到同一类，通过混淆矩阵计算准确率、召回率和ARI。公式中TP(真阳性)和FP(假阳性)分别表示真正属于同一物种的contigs聚在同一类和不同类中的数目；FN(假阴性)和TN(真阴性)分别代表不同物种的contigs聚集成同一类和不同类的数目。
下面通过具体实施例对本发明一种基于参考物种标签约束的宏基因组序列深度聚类方法做进一步的解释说明：
实施例1
本实施例中，测试集包括Sharon模拟数据集和Strain模拟数据集两个数据集，Sharon模拟数据集中包含101个物种的37628条contigs序列，是基于HMP项目的前96个数据集进行了模拟。Strain模拟数据集是测试不同物种分类分辨率的能力，包含20个物种的9401条contigs序列，这其中包括同一物种的菌株。20个物种分别为5个不同的大肠杆菌菌株、5种拟杆菌、5种来自不同梭菌属的菌株以及5种其他典型肠道细菌菌株组成。
利用本发明提供的基于物种标签约束的宏基因组深度聚类方法(Label-constrained deep clustering，LCDC)进行分析，同时将目前的同类分析方法COCACOLA、CONCOCT、MetaBAT和MaxBin2.0作为对照，分析结果如图5、图6所示。
101物种的模拟数据集是包含了96个样本的37628条重叠群序列，每个样本包含来自相同101个物种但相对频率不同的随机配对末端读数。由图5可知，COCACOLA在101物种的Sharon模拟数据集上表现的非常不稳定，ARI仅有62.3％；CONCOCT、MetaBAT和MaxBin2.0均超过了90％；LCDC在101物种数据集中ARI系数表现最优，达到97.2％。
对于Sharon数据集，使用人体肠道微生物库预训练的模型，此模型的结构网络设置为136*100*80*50*80*100*136的五个隐藏层的对称结构。图6为20个物种的模拟数据集中各方法的聚类结果，在20物种的Strain模拟数据集中，COCACOLA表现了最好的效果，ARI达到了92％；LCDC也依然表现不错，ARI为86.7％。在Strain模拟数据集上，64个样本中仅包含20个基因组，且5个基因组来自不同大肠杆菌菌株的同一个属。计算结果分别如表1、2所示，表3为构建预训练数据集的物种信息。
表1、101物种的Sharon模拟数据中个方法的聚类结果
| COCACOLA | 64.2 | 80.8 | 92.1 |
| CONCOCT | 94.1 | 94.9 | 97.7 |
| MaxBin2.0 | 95.8 | 96.8 | 98.2 |
| MetaBAT | 94.9 | 96.2 | 93.8 |
| LCDC | 97.2 | 98 | 97.2 |
表2、20物种的Strain模拟数据集中个方法的聚类结果
| COCACOLA | 92 | 95.3 | 91.8 |
| CONCOCT | 87.2 | 86.9 | 97.1 |
| MaxBin2.0 | 81.8 | 91.7 | 75 |
| MetaBAT | 88.3 | 95.3 | 85.2 |
| LCDC | 87 | 86.2 | 92.5 |
表3、构建预训练数据集的物种信息
综上所述，本发明提供了一种基于参考物种标签约束的宏基因组序列深度聚类方法(Label-constrained deep clustering，LCDC)，所述方法与目前的同类方法有明显不同，其设计了基于参考物种标签约束的深度学习预训练模型。由于所有的生态系统都存在微生物群落，而且不同环境下的微生物群落存在巨大的差异，相同环境下也会因为温度、酸碱度、光照、通气性等细微差距而使得微生物种类存在巨大差异。本发明利用微生物群落差异的特点，对不同环境下的微生物群落进行分析，建立了基于不同群落的已知物种的预训练数据库，构建预训练数据库时将每条4mer特征向量分为同一物种、相同属不同物种和不同属不同物种三种情况，并分别研究了三种情况下的样本间序列的4mer特征间的关系；建立了预训练模型的标签约束误差函数，并且LCDC使用群落已知标签的数据库进行预训练，针对不同的微生物群落构建不同预训练模型；在用户使用时，只需要针对不同的群落加载所需群落的预训练模型，重新加载模型仅仅等待几次微调步骤的迭代即可得到聚类结果。经过多个数据集上进行分析，并对比LCDC与其他多种方法的聚类结果，在真实数据集中，LCDC能够展现非常优秀的聚类性能。
应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。
Claims (10)
1.一种基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，包括步骤：
第一步，预训练步骤，包括：
1.1构建基于不同微生物群落的已知物种的预训练数据库；
1.2建立初始化模型；
1.3针对不同的微生物群落构建不同的预训练模型；
第二步，微调步骤，包括：
2.1计算待测微生物群落的数据集重叠群序列4mer频率，并归一化得到微调的输入特征频率Y；
2.2加载待测微生物群落的预训练模型以及参数；
2.3构建微调误差函数；
2.4确定聚类个数；
2.5微调模型；
2.6完成聚类，得到聚类结果，并根据聚类指标衡量聚类结果。
2.根据权利要求1所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述步骤1.1具体包括步骤：
a)下载不同微生物群落的已知物种的全基因组序列预训练数据集；
b)将每个物种的全基因组序列从随机起始位置截取随机长度的序列；
c)计算步骤b)中截取的每条序列的4mer频率特征，并进行归一化，得到不同微生物群落的宏基因组预训练4mer频率归一化特征X。
3.根据权利要求1所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述步骤1.2具体包括步骤：
a)建立具有对称结构的自编码器；
b)选取
函数作为激活函数，加入Dropout函数来调节模型参数和样本量之间的关系；
c)设置模型参数；
d)构建预训练误差函数，其计算公式为：
ERRORpre＝ERES+κELCN
其中，ERES表示重构误差，ELCN表示标签约束误差，κ表示用于平衡重构误差ERES和标签约束误差ELCN的超参数。
4.根据权利要求3所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述重构误差ERES的计算公式为：
其中，xi表示编码网络的输入，f(xi)表示编码网络的输出，g(f(xi))表示解码网络的输出，N1表示总样本的个数；
所述标签约束误差ELCN的计算公式为：
其中，Es表示衡量相同物种间的特征向量间的欧式距离，El表示衡量相同属不相同物种间的特征向量间的距离，Ed表示衡量不同属不同物种间的特征向量间的距离，n1、n2、n3为三种误差下累加的次数且满足
β、λ是标签约束相的超参数。
5.根据权利要求1所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述步骤1.3中，构建预训练模型具体包括步骤：
a)使用初始化后的网络模型以及参数；
b)加载测试集样本归一化特征X，并送入初始化后网络模型，计算重构误差和标签约束误差；
c)应用反向传播自适应矩估计方法，对不同变化的参数以自适应的学习率进行更新；
d)保存预训练的模型以及参数，定义为待测微生物群落的预训练模型。
6.根据权利要求1所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述步骤2.3具体包括步骤：
利用深度k-mean聚类的方法，将k-mean聚类的误差加入到预训练得到的待测微生物群落的预训练模型中，得到微调误差函数：
其中，ERES表示重构误差，ECLU表示聚类误差，N2表示待测微生物群落样本的个数，η表示用于平衡重构误差ERES和聚类误差ECLU的超参数。
7.根据权利要求6所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述重构误差ERES计算公式为：
其中，yi表示微调网络测试集的输入，g(f(yi))表示微调网络测试集的输出；
所述聚类误差ECLU计算公式为：
其中，x表示Ci类中的样本点，θi表示Ci类的聚类中心。
8.根据权利要求1所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述步骤2.4具体包括步骤：
a)定义Si变量为类内数据到簇质心的平均距离，代表簇Ci中各样本的分散程度，其计算公式为：
其中，|Ci|表示簇类Ci中样本点个数，x表示Ci类中的样本点，θi表示Ci类的聚类中心；
b)定义Mij为簇Ci与簇Cj的距离，其计算公式为：
Mij＝||θi-θj||2
其中，Mij表示簇Ci与簇Cj的质心θi和θj的距离；
c)应用DBI指数确定聚类个数k，其计算公式为：
其中，k为聚类个数变量，取值范围为[0.002*N2,0.005*N2]；
d)选择最优的聚类个数k0，其计算公式为：
其中，k为聚类个数变量。
9.根据权利要求1所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述步骤2.5具体包括步骤：
a)将测试样本送入网络，得到第一次隐藏表示；
b)将a)中得到的隐藏表示作为样本的特征使用k-mean聚类，聚类个数为k0，初始化中心的方法使用k-means++；
c)将b)中得到的聚类结果分批次随机取出，计算重构误差ERES和聚类误差ECLU，并根据重构误差ERES和聚类误差ECLU计算得到微调误差ERRORfine；
d)对每个样本的微调误差ERRORfine使用Adam优化器微调模型参数，得到微调后的模型。
10.根据权利要求1所述的基于参考物种标签约束的宏基因组序列深度聚类方法，其特征在于，所述步骤2.6具体包括步骤：
a)联合样本的覆盖率信息和网络中间层用于聚类，得到聚类结果
b)使用准确率、召回率和ARI来衡量聚类结果，其中，准确率的计算公式为：
召回率的计算公式为：
ARI的计算公式为：
其中，TP和FP分别表示真正属于同一物种的contigs聚在同一类和不同类中的数目，FN和TN分别表示属于不同物种的contigs聚集成同一类和不同类的数目。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111389111.2A CN114065866B (zh) | 2021-11-22 | 2021-11-22 | 一种基于参考物种标签约束的宏基因组序列深度聚类方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111389111.2A CN114065866B (zh) | 2021-11-22 | 2021-11-22 | 一种基于参考物种标签约束的宏基因组序列深度聚类方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114065866A true CN114065866A (zh) | 2022-02-18 |
| CN114065866B CN114065866B (zh) | 2024-04-30 |
Family
ID=80279297
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111389111.2A Active CN114065866B (zh) | 2021-11-22 | 2021-11-22 | 一种基于参考物种标签约束的宏基因组序列深度聚类方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114065866B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080113872A1 (en) * | 2004-04-29 | 2008-05-15 | Frank Henri Johan Schuren | Classification of Organisms Based on Genome Representing Arrays |
| CN106599618A (zh) * | 2016-12-23 | 2017-04-26 | 吉林大学 | 一种宏基因组重叠群的无监督分类方法 |
| US20180137243A1 (en) * | 2016-11-17 | 2018-05-17 | Resilient Biotics, Inc. | Therapeutic Methods Using Metagenomic Data From Microbial Communities |
| CN112466404A (zh) * | 2020-12-14 | 2021-03-09 | 浙江师范大学 | 一种宏基因组重叠群无监督聚类方法及系统 |
-
2021
- 2021-11-22 CN CN202111389111.2A patent/CN114065866B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080113872A1 (en) * | 2004-04-29 | 2008-05-15 | Frank Henri Johan Schuren | Classification of Organisms Based on Genome Representing Arrays |
| US20180137243A1 (en) * | 2016-11-17 | 2018-05-17 | Resilient Biotics, Inc. | Therapeutic Methods Using Metagenomic Data From Microbial Communities |
| CN106599618A (zh) * | 2016-12-23 | 2017-04-26 | 吉林大学 | 一种宏基因组重叠群的无监督分类方法 |
| CN112466404A (zh) * | 2020-12-14 | 2021-03-09 | 浙江师范大学 | 一种宏基因组重叠群无监督聚类方法及系统 |
Non-Patent Citations (3)
| Title |
|---|
| FAFU生信 小蘑菇: "聚类分析6—指示物种（数量生态学：R语言的应用 第四章） - 知乎", 《HTTPS://ZHUANLAN.ZHIHU.COM/P/371721658》, 12 May 2021 (2021-05-12), pages 1 - 17 * |
| TIZIAN SCHULZ等: "GraphTeams: a method for discovering spatial gene clusters in Hi-C sequencing data", 《BMC GENOMICS》, 31 December 2018 (2018-12-31), pages 60 - 95 * |
| 疯狂的AI学长: "基于不同物种建立不同的聚类模型, 使用时对模型进行微调", 《HTTPS://ZHUANLAN.ZHIHU.COM/P/268880550》, 27 October 2020 (2020-10-27), pages 1 - 7 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114065866B (zh) | 2024-04-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Segal et al. | Rich probabilistic models for gene expression | |
| Otwinowski et al. | Genotype to phenotype mapping and the fitness landscape of the E. coli lac promoter | |
| CN106202999A (zh) | 基于不同尺度tuple词频的微生物高通量测序数据分析协议 | |
| CN112215259A (zh) | 基因选择方法和装置 | |
| CN114065866B (zh) | 一种基于参考物种标签约束的宏基因组序列深度聚类方法 | |
| Bai et al. | KIMI: Knockoff inference for motif identification from molecular sequences with controlled false discovery rate | |
| Dudoit et al. | Loss-based estimation with cross-validation: Applications to microarray data analysis | |
| Sekula | OptCluster: an R package for determining the optimal clustering algorithm and optimal number of clusters. | |
| Qiao et al. | Poisson hurdle model-based method for clustering microbiome features | |
| CN115579068A (zh) | 一种基于预训练和深度聚类的宏基因组物种重建方法 | |
| CN112967755A (zh) | 一种面向单细胞rna测序数据的细胞类型识别方法 | |
| CN108182347B (zh) | 一种大规模跨平台基因表达数据分类方法 | |
| CN111755074A (zh) | 一种酿酒酵母菌中dna复制起点的预测方法 | |
| CN113113137B (zh) | 基于最大相关最小冗余和改进花授粉算法的特征选择方法 | |
| Ramchandran et al. | Cross-cluster weighted forests | |
| Deng | Algorithms for reconstruction of gene regulatory networks from high-throughput gene expression data | |
| CN116364195B (zh) | 一种基于预训练模型的微生物遗传序列表型预测方法 | |
| Wu et al. | Cluster analysis of dynamic parameters of gene expression | |
| Walker | Iterative Random Forest Based High Performance Computing Methods Applied to Biological Systems and Human Health | |
| Marchetti-Bowick | Structured Sparse Regression Methods for Learning from High-Dimensional Genomic Data | |
| XING et al. | HANDLING HIGHLY CORRELATED GENES OF SINGLE-CELL RNA SEQUENCING DATA IN PREDICTION MODELS | |
| Essinger et al. | Neural network-based taxonomic clustering for metagenomics | |
| Hussami et al. | A component lasso | |
| Blein-Nicolas et al. | Nonlinear network-based quantitative trait prediction from biological data | |
| Banerjee et al. | Classification and clustering analysis of pyruvate dehydrogenase enzyme based on their physicochemical properties |
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 |