CN115579068A

CN115579068A - 一种基于预训练和深度聚类的宏基因组物种重建方法

Info

Publication number: CN115579068A
Application number: CN202211069609.5A
Authority: CN
Inventors: 刘富; 宋文智; 刘云; 苗岩
Original assignee: Jilin University
Current assignee: Jilin University
Priority date: 2022-08-31
Filing date: 2022-08-31
Publication date: 2023-01-06

Abstract

本发明提供了一种基于预训练和深度聚类的宏基因组物种重建方法。基于预训练和深度聚类的宏基因组物种重建方法，设计了基于图卷积神经网络联合Focal Loss损失函数的词嵌入特征提取模型以及基于LSTM自编码器联合改进的FCM算法的深度聚类模型。本发明构建的深度聚类算法模型，将深度学习与聚类结合在一起，重构误差与聚类误差同步优化，进一步提升二者性能，计算量也较小。在用户使用时，只需要针对所选取的数据集的大小及序列长度对整个模型的参数进行调整，重新运行模型即可得到聚类结果，大大提高了准确度与便利性，能够得到更加优秀的聚类结果。相较于其他方法，本发明所发现的未知物种完整度更高，污染度更低。

Description

一种基于预训练和深度聚类的宏基因组物种重建方法

技术领域

本发明涉及生物信息学分析领域，尤其涉及一种基于预训练和深度聚类的宏基因组物种重建方法。

背景技术

微生物个体微小，包括细菌、真菌及一些小型的原生生物等，是地球上种类最大、数量最多、分布最广的生物群。微生物之间的相互作用对人类健康和疾病等方面至关重要。长期以来，对微生物的研究都主要围绕着可培养的单一物种为对象而获得的研究成果。事实上，环境中可培养的微生物占自然界微生物总数的极少部分，经研究证实，1g土壤中存在10⁷个细菌物种，其中能被培养的仅占1％～10％左右，大多数的微生物种类与功能至今还没有被发现。

如今，下一代测序方法(Next-generation sequencing，NGS)，也就是高通量测序方法，能够详细分析环境样本中所有微生物的基因组，而不仅仅针对那些提前培养的微生物展开研究。宏基因组学则利用新一代测序技术在不经过实验室培养的情况下，能够获取环境中绝大部分的遗传物质。微生物群落聚类分析则是宏基因组研究中最具挑战性的任务之一。与传统的测序方法不同，第二代测序技术所得到的原始数据一般为大量的、长度较短的、物种间相对丰度不均匀的、来源于多种微生物的DNA片段。由于环境中众多物种的宏基因组DNA片段是相互交织的，故将其高效聚类和精确装箱是分析每个物种详细信息和对应功能的关键步骤。而且，每个DNA片段都包含较少量的碱基和特征，序列的复杂性和物种间相对不平衡性增加了后续聚类和拼接的难度。DNA片段的复杂性严重阻碍了宏基因组序列的准确研究和高效应用，增加了聚类模型构建和优化的难度。因此，如何高效、准确地聚类宏基因组DNA片段是目前需要解决的关键问题。大多数研究方法都是先将原始数据DNA片段进行组装，围绕组装后所得到的重叠群(contigs)展开分析的，而重叠群(contigs)具有长度较长、特征数量较多等特性，其数据集复杂性与不平衡性比DNA片段低。目前，虽然已有研究人员直接针对从环境中提取的宏基因组DNA片段的聚类和拼接算法提出了些许方法，但其性能还有待进一步提高。

因此，现有技术还有待改进。

发明内容

鉴于上述现有技术的不足，本发明的目的在于提供一种基于预训练和深度聚类的宏基因组物种重建方法，旨在解决因宏基因组DNA片段数据长度较短、特征较少的复杂性以及物种间相对丰度不均匀等特性带来的宏基因组DNA片段聚类与拼接难的问题。

本发明的技术方案如下：

本发明提供一种基于预训练和深度聚类的宏基因组物种重建方法，其中，包括步骤：

第一步，特征提取步骤，包括：

提取不同环境微生物的宏基因组原始数据集，并对原始数据集进行预处理，所述原始数据集包括不同物种的DNA序列特征；

针对预处理后的数据集构建词嵌入模型；

构建模型误差函数；

利用构建的模型误差函数对词嵌入模型进行训练并调节参数；

保存输出的特征向量矩阵；

第二步，深度聚类步骤，包括：

构建深度LSTM自编码器联合聚类模型，所述LSTM自编码器包括编码器和解码器两部分，编码器对输入的特征向量时间序列数据进行特征学习，解码器利用当前隐藏层状态及网络参数对数据特征进行重构；

通过LSTM自编码器联合FCM算法聚类模型来构建模型聚类损失函数；

输入待测试的微生物的宏基因组数据集，利用模型重构损失函数及聚类损失函数计算模型的整体损失误差；

计算并分析模型的聚类性能指标；

调整模型参数，得到模型最优的聚类性能；

对待测试的数据集中的未知物种进行聚类、组装，并利用软件分析组装后的重叠群的完整度与污染度，输出未知物种的聚类结果。

所述的基于预训练和深度聚类的宏基因组物种重建方法，其中，所述对原始数据集进行预处理具体包括步骤：

a)下载微生物群落的宏基因组序列数据集；

b)根据数据集中存储的每条DNA片段的质量值信息，通过质控软件工具进行优化，过滤掉低质量序列；

c)将步骤b)过滤后剩余的序列中的碱基N替换成A、G、C、T中的一种；

d)利用BLAST工具比对步骤c)处理后的宏基因组全部序列的标签，根据引物和索引序列区分不同样本；

e)根据BLAST比对结果，过滤掉不属于同一物种的DNA片段对，将过滤后的序列作为词嵌入模型的输入。

所述的基于预训练和深度聚类的宏基因组物种重建方法，其中，所述低质量序列包括序列中连续存在超过5个碱基N的序列。

所述的基于预训练和深度聚类的宏基因组物种重建方法，其中，所述构建词嵌入模型具体包括步骤：

a)利用滑窗法对预处理后的DNA片段序列进行切割，将序列转化为重叠的固定长度k-mer序列；

b)将切割后的序列数据库构建成拓扑图；

c)通过k-mer序列的共现信息来建立两个k-mer节点之间的边；

d)将构建好的拓扑图输入至两层图卷积神经网络中，构建词嵌入模型。

所述的基于预训练和深度聚类的宏基因组物种重建方法，其中，所述步骤构建模型误差函数中，采用Focal Loss作为特征提取模型的误差函数。

所述的基于预训练和深度聚类的宏基因组物种重建方法，其中，所述编码器和解码器分别包括两层LSTM和两层全连接层，编码器和解码器的中间构建了一层聚类层。

所述的基于预训练和深度聚类的宏基因组物种重建方法，其中，所述LSTM的神经元由若干个递归连接的记忆区块构成，每个记忆区块包括1个记忆单元以及3个逻辑单元，所述逻辑单元包括遗忘门、输入门和输出门。

所述的基于预训练和深度聚类的宏基因组物种重建方法，其中，所述步骤计算并分析模型的聚类性能指标中，所述聚类性能指标包括精确率、召回率以及调整兰德系数。

一种存储介质，其中，所述存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如如上任意所述的基于预训练和深度聚类的宏基因组物种重建方法中的步骤。

一种终端设备，其中，包括处理器，适于实现各指令；以及存储介质，适于存储多条指令，所述指令适于由处理器加载并执行如上任意所述的基于预训练和深度聚类的宏基因组物种重建方法中的步骤。

有益效果：本发明提供了一种基于预训练和深度聚类的宏基因组物种重建方法。基于宏基因组DNA片段的深度LSTM自编码器联合聚类方法，设计了基于图卷积神经网络联合Focal Loss损失函数的词嵌入特征提取模型以及基于LSTM自编码器联合改进的FCM算法的深度聚类模型。本发明将所选取的宏基因组数据集中每条DNA片段用滑窗法以k-mer的形式编码，应用两层GCN网络作为词嵌入模型训练高维特征，将宏基因组中没有比对标签的片段用训练后得到的词嵌入结果表示，作为后续深度聚类的数据样本。本发明构建了一种深度LSTM自编码器联合聚类算法模型，对上述处理后的数据集进行聚类分箱，将深度学习与聚类结合在一起，重构误差与聚类误差同步优化，相比于其他算法，可以进一步提升二者性能，同时计算量也较小。最后，利用拼接软件对聚类后的片段按照类别进行组装，并用check-m工具计算组装后的重叠群(contigs)的完整度与污染度。在用户使用时，只需要针对所选取的数据集的大小及序列长度对整个模型的参数进行调整，重新运行模型即可得到聚类结果，大大提高了准确度与便利性。经过真实数据集的验证以及与现有的其他多种方法的聚类结果进行比较，结果证实，本发明能够得到更加优秀的聚类结果。不仅如此，在发现未知物种的层面上，本发明所发现的未知物种相较于其他方法而言，完整度更高，污染度更低。

附图说明

图1为本发明实施例中一种基于预训练和深度聚类的宏基因组物种重建方法的整体模型流程图。

图2为本发明实施例中DNA-FLGCN模型结构示意图。

图3为本发明实施例中滑窗法构建词嵌入特征向量矩阵流程图。

图4为本发明实施例中LSTM网络神经元结构示意图。

图5为本发明实施例中LSTM自编码器解码器的处理方式示意图(语义向量仅作初始状态参与运算)。

图6为本发明实施例中LSTM自编码器解码器的处理方式示意图(语义向量参与所有时刻运算)。

图7为本发明实施例中深度LSTM自编码器联合聚类模型结构示意图。

图8为本发明实施例中不同方法下SRR492190测试集聚类结果的比较。

图9为本发明实施例中不同方法下SRR492190未知物种片段聚类拼接后重叠群的完整度比较。

图10为本发明实施例中不同方法下SRR492190未知物种片段聚类拼接后重叠群的污染度比较。

具体实施方式

本发明提供一种基于预训练和深度聚类的宏基因组物种重建方法，为使本发明的目的、技术方案及效果更加清楚、明确，以下对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明实施例提供一种基于预训练和深度聚类的宏基因组物种重建方法，包括步骤：

S100、特征提取步骤；

S200、深度聚类步骤。

由于宏基因组有很多序列与NCBI官方核酸序列数据库中的序列相似度较低，导致无法比对出这些序列的标签，至今为止大多数方法都是基于有比对标签结果的序列片段而展开的，未对没有比对标签结果的序列进行详细分析。本发明实施例基于宏基因组DNA片段的深度LSTM自编码器联合聚类方法与目前的同类方法存在明显不同，设计了基于图卷积神经网络联合Focal Loss损失函数的词嵌入特征提取模型以及基于LSTM自编码器联合改进的FCM算法的深度聚类模型。图1所示为本发明实施例所述基于预训练和深度聚类的宏基因组物种重建方法的整体模型流程图。本发明实施例直接提取环境中通过高通量测序技术所得到的宏基因组原始数据集(reads)，对于此类特征较少、长度较短、物种间相对丰度极不均匀的数据集展开详细研究。利用构建好的词嵌入模型提取数据集特征，得到特征向量矩阵，构建边训练边聚类的深度自编码器联合聚类模型对数据集聚类并进行详细分析，在提升聚类性能的同时，有效降低模型的时间复杂度。本发明实施例主要针对所选取的宏基因组数据集中未知标签的序列进行分析，将其用训练词嵌入模型后得到的特征向量表示，输入至构建好的深度聚类模型中并进行调参，从而得到较优的聚类性能指标，并按照聚类后的结果进行拼接，分析拼接后的重叠群(contigs)的完整度与污染度。

在一些实施方式中，S100具体包括以下步骤：

S101、提取不同环境微生物的宏基因组原始数据集，并对原始数据集进行预处理，所述原始数据集包括不同物种的DNA序列特征；

S102、针对预处理后的数据集构建词嵌入模型；

S103、构建模型误差函数；

S104、利用构建的模型误差函数对词嵌入模型进行训练并调节参数；

S105、保存输出的特征向量矩阵。

在一些实施方式中，所述步骤S101中，对原始数据集进行预处理具体包括步骤：

a)下载微生物群落的宏基因组序列数据集；

宏基因组中存在大量不同物种的DNA序列，因此提取每条DNA序列的特征并按照其物种进行分类是宏基因组片段序列分析的重要一步。DNA的序列特征是指A、T、G、C四种核苷酸的排列组合顺序，如果两个DNA序列具有相似的组合，则这两个DNA序列的相似性程度较高，其来源于同一物种的可能性就较高。

具体的，本发明实施例通过NCBI网站下载不同微生物群落的宏基因组序列数据集，例如人体肠道微生物宏基因组序列数据集。作为一个具体的实施例，以SRR492190数据集为例。高通量测序技术测序得到的双端序列数据，以fastq或fasta格式存储。以fastq格式为例，下载后的高通量测序序列分为fq1和fq2两个文件，其中fq1文件里的序列被称为read1，fq2文件里的序列被称为read2，相对应的read1和read2序列为两条反向互补链，每条片段的长度均为100bp。

具体的，本发明实施例根据文件中存储的每条DNA片段(reads)的质量值信息通过质控软件fastp工具进行优化。利用fastp等数据质控软件过滤低质量的序列(例如连续大量N碱基存在的序列)。更为具体的，本发明实施例规定，当一条序列中连续存在超过5个碱基N时，将该序列也被视为低质量序列并将其过滤。

具体的，将过滤后剩余的序列中的碱基N随机替换成A、G、C、T任意一种。若下载的高通量测序文件中序列存在碱基N，则说明该序列碱基N所在的位置没有得到明确的测序结果，故而本发明实施例将其随机替换成A、G、C、T任意一种。

具体的，利用BLAST工具比对得到处理后的宏基因组全部序列的标签，根据引物和索引序列区分不同样本。由于fq1与fq2文件中DNA片段为双端测序结果，故每一对反向互补的DNA片段应属于同一物种。检测BLAST比对结果，过滤掉不属于同一物种的DNA片段对，将过滤后的fq1文件中的序列作为词嵌入模型的输入。

在一些实施方式中，所述步骤S102中，构建词嵌入模型具体包括步骤：

b)将切割后的序列数据库构建成拓扑图；

c)通过k-mer序列的共现信息来建立两个k-mer节点之间的边；

本发明实施例建立了一种基于DNA片段(reads)的词嵌入模型(DNA-FLGCN)，在Text GCN的基础上进行改进与创新。

具体的，利用滑窗法对预处理后的fq1文件中的序列进行切割，将序列转化为重叠的固定长度k-mer序列。假设每条序列长度为n，则该序列中共有(n-k+1)个窗口。例如，将ATGCCGTAT转化为5个5-mer序列:{ATGCC,TGCCG,GCCGT,CCGTA,CGTAT}。

本发明实施例所提出的DNA-FLGCN模型结构图如图2所示。

具体的，将切割后的序列数据库构建成拓扑图，该图的节点由整个文档和所有k-mer序列构成，显式地对全局的单词共现进行建模，便于后续用图卷积神经网络(GCN)进行训练。

具体的，通过k-mer序列的共现信息来建立两个k-mer节点之间的边。每个k-mer节点和文档节点之间的边是通过该k-mer的频率和k-mer所在的文档频率来建立的。两个k-mer节点之间的边的权重定义如下：

其中TF-IDF(i,j)表示文档节点和每个k-mer序列节点之间的边的权重，PMI(i,j)为两个k-mer序列节点间的权重。当PMI(i,j)为正数时，表示序列i和序列j有较强的语义关联性，当PMI(i,j)为负数的时候，则表示词汇i，j语义关联性较低。

具体的，将构建好的拓扑图输入至两层图卷积神经网络(GCN)中，构建词嵌入模型。第一层图卷积神经网络(GCN)输出大小为词嵌入向量训练后的维度，第二层图卷积神经网络(GCN)输出即为标签的总类别个数。在所构建的词嵌入模型中，测试精度随着窗口尺寸的增大而增大，但过小的窗口并不能充分产生全局词共现信息，窗口尺寸过大可能使得联系度较低的k-mer序列之间添加边；同时，词嵌入维度大小也为影响整个模型的重要因素，低维的词嵌入维度不能较好地将完整的标签信息传播到所构建的图中，词嵌入维度过高可能会降低分类性能，并且会消耗较多的训练时间。优选的，本发明实施例设置窗口尺寸为10，词嵌入维度为10。模型选取

函数作为激活函数，加入Dropout函数防止模型过拟合，通过该函数调节模型参数和样本数量之间的关系，设置Dropout函数的参数为0.5。

在一些实施方式中，所述步骤S103中，采用Focal Loss作为分类模型的误差函数。

现实生活中大多数宏基因组数据集极度不平衡，物种间相对丰度差异较大，仅使用传统的交叉熵作为损失函数并不能将其进行准确分类，故需进一步改进损失函数。根据样本分辨的难易程度给对应的损失添加权重。其中，分类置信度接近1或接近0的样本称为易分辨样本，其余的称之为难分辨样本。考虑多方面因素，本发明实施例采用Focal Loss作为分类模型的误差函数。

具体的，Focal Loss误差函数的计算公式为：

FL(p_t)＝-α_t(1-p_t)^γlog(p_t) (2)

式中，p_t为当前样本所属的类别在总数据集中所占的概率分布，在多分类问题中，即为将图卷积神经网络训练后得到的预测标签通过softmax层后输出的概率。α_t与γ是两个超参数，α_t主要作用是解决正负样本的不平衡问题(α_t∈[0,1])，样本所在类别数量越多，α_t越小，目的在于降低高丰度物种训练损失，使后续分类与优化的注意力集中在低丰度及极低丰度物种上。γ则是解决难易样本的不平衡问题，其能够有效降低负样本的损失，当γ越大时，概率越大的简单样本的Focal Loss值可以降的越低。优选的，γ取值为2时效果最佳。

本发明实施例提出的DNA-FLGCN方法对每个k-mer序列进行编码，即使用一个特征向量表征一个k-mer，特征向量的每个元素都是对该k-mer某一特征的量化描述。k-mer中k值的选择会影响后续基因组聚类及组装。较小的k-mer将减少所构建图的边，增加所有k-mer重叠的机会，减少存储DNA片段所需内存，但会面临多顶点通向单个k-mer序列的风险，会丢失数据集中一部分信息；较大的k值有助于基因组的构建，但会导致DNA片段的脱节，在组装过程中会生成大量长度较短的重叠群(contigs)。

优选的，本发明实施例取k值为13，如图3所示，利用滑窗法构建词嵌入特征向量矩阵，其中L为DNA片段序列长度，以本发明实施例所用数据集SRR492190为例，L值为100，每条序列可通过滑窗法分割成88个13-mer序列。特征向量的长度依具体情况而定，特征元素越多则对该单词表征得越全面。优选的，设置特征向量长度为10，序列中每个k-mer序列都可以使用对应的88×10的向量特征来表示。该方法所生成的词向量具有很好的语义特性，并且其是神经网络在自然语言处理领域应用的结果，利用深度学习方法来获取词语的分布表示，可以用于文本分类、情感计算、词语构建等自然语言处理任务。这种特征表征的优势在于可以清楚表明不同k-mer序列之间的相似程度，且处理后的序列特征数量明显增多，更加有利于后续自编码器深度学习与聚类。

在一些实施方式中，所述步骤S104中，利用构建的模型误差函数对词嵌入模型进行训练并调节参数具体包括：将预处理后的fq1文件中的数据集作为词嵌入模型的输入，其中随机选取75％的序列作为训练集，其余25％的序列作为测试集，将其输入至上述构建好的词嵌入模型中。在训练词嵌入模型时，由于需要详细地构建每条序列信息的相关性，故须在一次迭代中将训练集全部数据一次性输入到模型中，并通过测试集计算每一次迭代后词嵌入训练模型的精确率(Precision)、召回率(Recall)、调整兰德系数(Adjusted RandIndex，ARI)等性能指标。以本发明实施例所用数据集SRR492190为例，本发明优选设置自适应矩估计(Adaptive moment estimation，Adam)优化器的学习率lr_rate＝0.15，训练迭代周期epoch＝600。对于不同的数据集，需要通过调整自适应学习率及训练迭代周期得到更佳的性能指标。

在一些实施方式中，所述步骤S105中，保存上述词嵌入模型输出的特征向量矩阵，为下一步宏基因组DNA片段深度聚类详细分析作基础。

在一些实施方式中，S200具体包括以下步骤：

S201、构建深度LSTM自编码器联合聚类模型，所述LSTM自编码器包括编码器和解码器两部分，编码器对输入的特征向量时间序列数据进行特征学习，解码器利用当前隐藏层状态及网络参数对数据特征进行重构；

S202、通过LSTM自编码器联合FCM算法聚类模型来构建模型聚类损失函数；

S203、输入待测试的微生物的宏基因组数据集，利用模型重构损失函数及聚类损失函数计算模型的整体损失误差；

S204、计算并分析模型的聚类性能指标；

S205、调整模型参数，得到模型最优的聚类性能；

S206、对待测试的数据集中的未知物种进行聚类、组装，并利用软件分析组装后的重叠群的完整度与污染度，输出未知物种的聚类结果。

在一些实施方式中，所述步骤S201中，所述编码器和解码器分别包括两层LSTM和两层全连接层，编码器和解码器的中间构建了一层聚类层。

在一些具体的实施方式中，LSTM神经元由若干个递归连接的记忆区块构成。具体的，每个记忆区块包括记忆单元以及3个逻辑单元所述，逻辑单元包括遗忘门、输入门和输出门。

长短期记忆网络(LSTM)神经元结构如图4所示，其是由一系列递归连接的记忆区块构成，每个记忆区块包括了记忆单元以及3个逻辑单元，通过引入这3个逻辑单元控制信息的传递，分别命名为：遗忘门、输入门和输出门。图中每个浅色圆形为对向量做出相应操作，每个矩形代表一个神经网络层，其内部字符代表相应神经网络所使用的激活函数。图中从单元中贯穿的线为神经元的隐藏状态，σ代表sigmoid函数，即

当x趋近于负无穷时，f(x)趋近于0，是作用于神经元的非线性函数。c_t表示神经元在经过t时刻后的记忆，其涵盖了在t+1时刻前神经网络对其所有输入信息的“总结概括”，c_t-1表示神经元在前一时刻的记忆。记忆单元表示神经元状态的记忆，输入门和输出门分别用来接收和输出参数，输入门用tanh函数提取有效信息，用

表示，使用sigmoid函数来控制哪些记忆要放至单元状态，用i_t表示，为每个分量做出评级，评级越高的则会有更多的记忆输入至LSTM记忆单元中。输出门将当前输入值与前一时刻输出值整合后的向量用sigmoid函数提取信息，再将当前记忆单元通过tanh函数压缩至区间(-1,1)中，再将处理后的单元状态与sigmoid函数处理后的结果点乘即可得到LSTM神经网络在t时刻的输出；遗忘门则用来控制是否保留当前隐藏层节点存储的历史信息，LSTM神经网络会依据新的输入及前一时刻输出决定遗忘掉之前的哪部分记忆，其中将新的输入x_t和前一步的输出h_t-1整合为单独向量，通过sigmoid层点对点的乘在当前单元状态上，用o_t表示，sigmoid函数将输入压缩到(0,1)区间中，故若整合后的向量中某个分量在通过sigmoid层后变至0，则相应单元状态在对位相乘后对应分量亦会变成0，即“遗忘”该分量上的信息，反之则表示保持完整记忆，h_t为当前时刻的输出。长短期记忆网络(LSTM)可长期记忆重要信息，并可随着输入进行动态调整，其中f_t为sigmoid神经层的输出向量。长短期记忆网络(LSTM)有了读取、重置和更新历史信息的能力。为了最小化训练误差，利用梯度下降法可修改每次训练时的网络权重。

在一些实施方式中，所述步骤S201中，LSTM自编码器包括编码器和解码器两部分，编码器对输入的特征向量时间序列数据进行特征学习，解码器利用当前隐藏层状态及网络参数对数据特征进行重构。

Seq2seq模型属于LSTM自编码器结构的一种，主要思想是利用两个循环神经网络(RNN)，其中一个作为编码器，另一个作为解码器。编码器负责将输入序列进行降维，输出指定长度的向量，解码器则根据语义向量生成指定的序列。解码器的处理方式有两种，如图5和图6所示，其中h₀,h₁,h₂,h₃为编码器的状态向量，h₀’,h₁’,h₂’,h₃’为解码器的状态向量，x₁,x₂,x₃,x₄为输入向量，y₁,y₂,y₃为输出向量。由图可知，其中一种是语义向量C只作为初始状态参与运算，解码器后面部分的运算都与C无关；另一种处理方式是语义向量C参与序列所有时刻的运算。

LSTM自编码器的优势在于其解决DNA片段在训练过程中所存在的梯度消失和梯度爆炸的问题，并且能在不影响模型训练的前提下，可以有效地处理不定长序列。该模型可指定每一个批次中各个样本的有效序列长度，这样一条序列在有效长度内其状态值和输出值原理不变，但超过有效长度的部分的状态值将不会发生改变，而输出值都将变为零向量。输出值和状态值的计算不仅依赖当前时刻的输入值，也依赖于上一时刻的状态值。其内部原理是利用一个mask matrix矩阵标记有效部分和无效部分，无效部分不会造成反向传播时对参数的更新，同样也不会影响后续误差函数的计算。

假设输入数据序列为X_i，X_i＝{x_i1,x_i2,L,x_in}，对于每个序列X_i，此时t∈{1,2,L,n}所对应的第i列的隐藏层状态向量为：

式中

为第i个编码单元在(t-1)时刻的输出状态向量，

为输入向量，W、R为m×d、m×m的稀疏权重矩阵；函数k(·)通常设置为tanh激活函数，将X_i中的每个列向量输入到编码器部分，输出结果为：

式中

为第i个编码单元在t时刻的输出，α为编码器部分的参数，

同样为tanh激活函数。在高维数据的训练过程中，池化(pooling)是必要的，主要用于降维。池化方法包括随机池化、均值池化和最大池化等。本文使用最大池化方法进行训练，相比于其他方法，其更好地保留了数据原始特征，降低模型计算复杂度。经过池化后，将其输入到解码器部分，重构结果为：

其中，

为重构数据，

为第i个编码单元在t时刻的输出，

为解码器部分隐藏状态向量，

和ρ(·)通常设置为tanh函数，从而由公式

计算得到重构误差，其中，

为重构数据，

为原始数据。

在一些实施方式中，本发明实施例设定构建的深度聚类模型中，编码器和解码器分别包括两层长短期记忆网络(LSTM)和两层全连接层，在编码器和解码器的中间构建一层聚类层，实现边学习边聚类功能，在提高模型聚类性能的同时，降低时间复杂度，并且根据不同数据集适当更改模型的层数及参数。

在一些实施方式中，所述步骤S202中，通过LSTM自编码器联合FCM算法聚类模型来构建所述模型聚类损失函数。

现阶段对数据集进行特征学习与聚类的研究中，大多数方法都是在对其聚类或分箱之前对序列特征进行深度学习，重构损失与聚类损失交叉优化，但是聚类效果并不是十分稳定，并且模型的时间复杂度较高。利用无监督聚类算法联合优化深度神经网络已经成为一个活跃的研究领域，结合深度学习的探索方法模型包括基于区分和软聚类、K-means聚类、子空间聚类、Student’s-t分布、图聚类、基于Kullback-Leibler(KL)发散聚类等等。对于上文所述的宏基因组短片段聚类难的问题，本发明实施例提出了一种LSTM自编码器联合改进的FCM算法聚类模型，将深度学习与聚类结合在一起，重构误差与聚类误差同步优化，可以进一步提升二者性能，计算量也较小。

本发明实施例通过在自编码器网络产生的特征空间中表示数据，并提取隐藏层输出向量，利用一个聚类层将数据特征从高维空间映射到低维空间，迭代优化聚类目标并计算Student’s-t分布，以测量隐藏层与聚类中心间的相似度。模型中N个输入数据表示为X＝{x₁,…,x_i,…,x_N},x_i∈R^d，自编码器隐藏层输出定义为Z＝{z₁,…,z_i,…,z_N},z_i∈R^c，其中d、c分别为输入层和隐藏特征层的维数。将z_i输入到所构建的聚类层，实现基于Student’s-t分布的聚类，公式如下所示：

具体来说，上式定义了Student’s-t分布后的聚类层的输出概率分布q_ij，其中μ＝[μ₁,…,μ_j,…,μ_c]为训练后的聚类层权重，c为类别数，j为聚类层中的j_th神经元。q_ij为x_i对第j个簇的隶属度，μ为聚类中心或代表Student’s-t分布，α为超参数，通常设为1。

传统的模糊C-means聚类算法处理不平衡数据集时会存在“均匀效应”，聚类过程中会倾向于产生尺寸大小相同的簇，使得聚类效果较差，精确率与召回率等性能指标较低。假定数据集仅包括两个类，则FCM目标函数可调整为：

其中，n表示样本总数，n₁和n₂分别代表每个类所含的样本数，v₁和v₂代表两个类的聚类中心。由于输入数据x_i，x_j为常量，故

为定值，最小化L_FCM等价于最大化2n₁n₂||v₁-v₂||²。若||v₁-v₂||²是定值，则当且仅当n₁＝n₂＝n/2时，n₁n₂取最大值，即L_FCM取得最优解。此时，模糊C-means产生均匀效应。

结合上述问题，本发明实施例将深度LSTM自编码器与改进后的FCS算法相结合，模型结构图如图7所示。令Z＝{z₁,…,z_i,…,z_N},z_i∈R^c为自编码器隐藏层的向量表示，其中c为隐藏层的特征维度，N为输入的数据量。将得到的隐藏层输出向量通过自己所构建的模糊聚类层，以Z为输入，提供对应的模糊隶属度矩阵q_ij：

式中m为模糊系数，一般取值为2，β为超参数，用于平衡类内距离与类间距离，μ＝[μ₁,…,μ_j,…,μ_c]为训练后的聚类层权重，c为类别数，j为聚类层中的第j个神经元。

定义为隐藏层向量的平均值，即

z_i为数据低维特征向量。与传统FCS算法相似，隶属度公式由损失函数推导而得，整个深度聚类模型的聚类损失函数为：

其中，将η_j如公式(9)所示，统一为

保证隐藏层输出向量在聚类时得到的模糊隶属度矩阵为正，m为模糊系数，

为模糊隶属度矩阵。

定义为隐藏层向量的平均值，即

将低维特征向量z_i的类间距离

定义为d₁，将类内距离||z_i-μ_j||²定义为d₂。一般情况下，在训练初期d₂的规模往往会过大。但随着训练的推进，d₂的规模逐渐减小，反而d₁的规模逐渐扩大，无法对d₂产生积极影响。若d₂-d₁＜0，则隶属度q_ij＜0，与约束目标不一致，且无法减缓均匀效应对聚类的影响。因此，为了保证d₁和d₂保持相同比例，d₁通过

归一化，其能够清晰反映d₁和d₂之间距离差异的敏感性。

基于Student’s-t分布的期望目标p_ij由下式定义，其不仅增加了q_ij的方差，同时将结果相对于第j个聚类中心进行归一化：

其中，隶属度q_ij与聚类层输出特征向量z_i相关，隶属度q_ij为构建的聚类层输出特征向量经过softmax函数后得到的概率分布结果。根据求得的p_ij和q_ij计算KL散度误差，通过计算KL散度使p_ij和q_ij更接近，并同时将损失函数最小化，KL散度误差函数如下式：

本发明实施例先利用K-means聚类对第一次训练迭代后得到的隐藏层向量初始化聚类中心，从第二次迭代开始，模型通过上述构建的聚类算法对聚类中心进行更新。针对簇个数的选取问题，常规的DBI等算法只适用于特征较多的重叠群(contigs)聚类问题，由于宏基因组DNA片段长度较短、特征较少，无法用这些算法确定聚类的簇具体个数，故在进行聚类时，需根据许所选取的宏基因组数据集中已知标签的序列物种分布大致确定未知物种序列的簇个数的范围。

本发明实施例采用Adam算法的优化器更新网络参数。Adam优化器相比于SGD算法及其他优化算法，具有更不容易陷于局部最优、速度更快，学习效果更为有效、尝试纠正其他优化技术中存在的一些问题，如学习率消失或是高方差的参数更新导致损失函数波动较大问题等优点。

在一些实施方式中，所述步骤S203中，将宏基因组已知标签的部分DNA片段作为测试集测试构建模型的聚类性能。

本发明实施例以SRR492190宏基因组数据集为例，将高丰度物种序列作为测试集聚类，该测试集共36778条，隶属于10个物种，最高丰度与最低丰度物种序列条数比为9:1，仍符合不平衡数据集。将测试数据集用滑窗法分割成若干个13-mer序列，通过上述词嵌入模型输出的特征向量矩阵表示，将其输入至构建的深度聚类模型。由于13-mer序列共有上千万种排列组合方式，可能会出现所构建的词嵌入特征向量矩阵不存在测试集中某条13-mer序列的情况，本发明实施例将此类13-mer序列用全零向量矩阵进行补充。以SRR492190测试集为例，每条DNA片段长度为100bp，用词嵌入特征向量表示后每条序列维度为88×10。本发明实施例所构建的深度聚类模型在聚类过程中实现边训练边聚类，故模型整体损失误差函数为重构误差与聚类误差损失之和，具体公式如下：

式中，L₁为类内距离指标，L₂为类间距离指标，L₃为重构误差函数，δ、γ为两个超参数(δ,γ∈(0,1))，根据聚类性能指标适当调整两个超参数的值，使得模型获得最优的聚类性能。

在一些实施方式中，所述步骤S204中，所述聚类性能指标包括精确率、召回率以及调整兰德系数。

对于传统的分类算法，通常使用分类的准确率对分类结果进行评价，准确率则为正确分类的样本个数除以样本总数。然而，对于不平衡数据而言，准确率并不适用与对其进行详细评价与分析。例如，一个包含两类的不平衡数据集，多数类的样本个数占总体的90％，如果仅将多数类的所有样本正确分类，即可得到90％的分类准确率。但该评价指标的结果是毫无意义的，因为其他10％的少数类序列并未被正确分类至相应的簇。基于此，本发明实施例利用精确率(Precision)、召回率(Recall)、调整兰德系数(Adjusted RandIndex，ARI)作为不平衡数据集聚类的评价性能指标，其皆由混淆矩阵计算而得。

混淆矩阵有效地反映实例类别划分与聚类结果的重叠程度，如表1所示。正类为少数类，负类为多数类，TP(True Positive)和FN(False Negative)分别代表正类样本被分到正类和负类的个数，FP(False Positive)和TN(True Negative)分别表示负类样本被分到正类和负类的个数。

表1混淆矩阵具体表现形式

具体的，

一般ε＝1。只有当召回率和精确率都较大时，F-value才能表示较好的聚类性能。因此，F-value可以评价聚类算法对于少数类的聚类性能。F-value定义为：

具体的，通过调整兰德系数(Adjusted Rand Index，ARI)衡量重叠程度，定义为：

在多分类评价指标中，考虑到存在样本的序列长短不一致性、分箱个数与物种真实标签种类个数不同、真实物种间的不平衡性等因素，本发明实施例对多标签聚类结果进行多角度分析，具体步骤包括：

a)对于序列长度不一致的情况，计算所述三个性能指标作为最终多标签聚类结果的总体评价指标。优选的，调整兰德系数(Adjusted Rand Index，ARI)用于聚类算法的比较时，取值范围为[-1,1]，精确率(Precision)和召回率(Recall)性能指标的取值范围为[0,1]，越接近1表示聚类效果越好。

b)分析物种间聚类的情况，对同一个属的多个物种进行聚类时，由于物种间存在相似的序列信息，难以将它们准确聚类，因此使用精确率和召回率分析聚类性能，更能表现出聚类效果的好坏。在属内分析得到更高的精确率及召回率，即说明此深度聚类算法模型对属内的不同物种间更敏感。

c)对于不平衡性较为明显的物种间序列的聚类性能时，结合上述所有性能指标统一分析，共同衡量聚类效果，避免分析不全面导致不准确的评价。

在一些实施方式中，所述步骤S205中，将所构建的测试集输入至深度LSTM自编码器联合聚类模型，根据多次迭代训练及聚类得到的性能指标合理调整模型参数。

在一些具体的实施方式中，针对SRR492190宏基因组中测试数据集，将深度聚类模型参数设定如下：编码器和解码器各有两层长短期记忆网络(LSTM)及两层全连接层，每条序列输入维度为88×10，两层长短期记忆网络(LSTM)输出维度分别为88×20及88×5，将每条序列经过长短期记忆网络(LSTM)后的输出维度合并为一维向量，随后通过两层全连接层网络，输出维度分别为200和64，将其通过上述构建的聚类层，使得每条序列的隐藏层输出维度与设定的簇类个数一致，便于在聚类过程中计算每条序列属于每个物种的概率，利用softmax函数得到概率最大值所在的物种，并将该序列聚类到该物种所在的类，编码器与解码器的结构为对称关系，解码器每一层网络的维度设置与编码器一致，使得最终自编码器的输出维度与输入维度一致，仍为88×10。模型整体损失误差函数中超参数δ＝0.5,γ＝0.9，迭代周期设置为500epoch，簇个数设置为10，每次迭代过程中，一个批次随机选取数据集中100条序列进行训练，经过一整个迭代周期后隐藏层特征向量矩阵更新，进行深度聚类并计算模型整体误差函数，当连续经过10个迭代周期后，模型的整体损失误差都没有下降时，输出聚类结果，计算并分析模型聚类性能指标。

在一些实施方式中，所述步骤S206中，对待测试的数据集中的未知物种进行聚类、组装，并利用软件分析组装后的重叠群的完整度与污染度，输出未知物种的聚类结果具体包括步骤：

a)将所述数据集的未知物种序列利用所构建的词嵌入模型输出的特征向量矩阵表示，并输入至构建的深度LSTM自编码器联合聚类模型中进行深度聚类；

b)将未知物种进行深度聚类后，利用Spades组装软件按照聚类结果进行组装；

c)利用check-m软件对组装后的重叠群进行分析，计算污染度与完整度。

本发明实施例整体的模型流程图如图1所示。将所选数据集的未知物种序列部分输入至所构建的深度LSTM自编码器联合聚类模型中。由于宏基因组数据集中，未知标签的序列与已知标签的序列的物种分布情况相似，而宏基因组DNA片段长度较短、特征较少，难以用常规算法确定类簇个数的大致范围，故需根据已知物种序列的分布确定未知标签序列的物种个数。

在一些具体的实施方式中，针对SRR492190宏基因组中测试数据集，该宏基因组已知标签序列共有12个高丰度物种(序列条数＞1000的物种)，序列总条数共占整个数据集约97％，其他物种视为极低丰度物种。故确定未知物种的簇类个数大致范围为[10，15]，最终本发明实施例根据该数据集将未知物种序列簇类个数设置为10。将未知物种序列通过上述词嵌入模型的特征向量矩阵表示，输入至深度LSTM自编码器联合聚类模型中，进行迭代调参并调整模型网络层数，最终当编码器与解码器各包括一层长短期记忆网络(LSTM)和三层全连接网络时，且自编码器结构设置为10*5*250*500*100*10*100*500*250*5*10的隐藏层的对称结构，学习率为0.001，迭代周期为200次，batch_size设为256，模型的整体聚类性能最佳。

将未知物种进行深度聚类后，利用Spades组装软件按照聚类结果进行组装，并利用check-m软件对组装后的重叠群(contigs)进行分析，计算其污染度与完整度，污染度越低且完整度越高说明该深度聚类模型效果更好。

本发明实施例还提供一种存储介质，所述存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如上任意所述的基于预训练和深度聚类的宏基因组物种重建方法中的步骤。

本发明实施例还提供一种终端设备，包括处理器，适于实现各指令；以及存储介质，适于存储多条指令，所述指令适于由处理器加载并执行如上任意所述的基于预训练和深度聚类的宏基因组物种重建方法中的步骤。

在一些实施方式中，所述存储介质中的逻辑指令可以通过软件功能单元的形式实现并作为独立的产品销售或使用。

存储介质作为一种计算机可读存储介质，可设置为存储软件程序、计算机可执行程序，如本公开实施例中的方法对应的程序指令或模块。处理器通过运行存储在存储介质中的软件程序、指令或模块，从而执行功能应用以及数据处理，即实现上述实施例中的方法。

存储介质可包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序；存储数据区可存储根据终端设备的使用所创建的数据等。此外，存储介质可以包括高速随机存取存储器，还可以包括非易失性存储器。例如，U盘、移动硬盘、只读存储(Read-Only Memory，ROM)、随机存取存储器(Random Access Memory，RAM)、磁碟或者光盘等多种可以存储程序代码的介质，也可以是暂态存储介质。

下面通过具体实施例对本发明一种基于预训练和深度聚类的宏基因组物种重建方法做进一步的解释说明：

实施例1

人体肠道微生物基因组研究计划(MetaHIT)是欧盟第七框架下属的子项目之一，目的是研究人类肠道微生物群落，进而了解功能和对人类健康的影响。本实施例中，选择了在人体肠道微生物基因组研究计划发表的文章中提取相关宏基因组数据集——Preborninfant gut metagenome-Carrol stool sample，以SRR492190数据集为例。

将整个数据集进行质控处理，并经过blast软件标签比对后共有1020842条DNA片段未得到对应标签，其中read1和read2各有510421条，本实施例主要目的为发现未知物种，且未知物种序列的分布应与已测得标签的物种序列分布一致。将得到标签的序列(共1523952条DNA片段，read1和read2各有761976条)过滤掉极低丰度物种，本实施例针对SRR492190数据集将序列条数小于1000条的物种定义为极低丰度物种，经过处理后剩余12个物种共1477442条DNA片段，其中最高丰度物种共有1206372条DNA片段，约占总数的81.65％。故该宏基因组数据集经过预处理后，仍具有极度不平衡性。

将本实施例与MBMC(2016)、Metaprob(2016)、MetaBin(2019)、Sparse Coding(2020)、Kexue Li,et al.(2020)、Metaprob2(2021)几种方法的聚类性能进行比较与分析。其中，Metaprob(2016)、Metaprob2(2021)解决了k-mer的可变分布以及不平衡的宏基因组片段聚类问题；MBMC(2016)利用马尔科夫链表示每个基因组，马尔可夫链描述了在状态空间中从一种状态过渡到另一种状态的随机过程；MetaBin(2019)是一种基于宏基因组不定长序列进行比较度量的方法；Kexue Li,et al.(2020)利用从数据集多个样本中获得的统计信息来减少欠聚类问题；Sparse Coding(2020)利用稀疏编码技术及弹性网络正则化实现潜在基因组恢复，并通过k-mer序列计数的稀疏性和非负性约束发现并恢复基因组的未知物种。本实施例相对于先前研究而言，将深度学习与聚类相结合，边训练边聚类的方案有效降低时间复杂度，将所选取的数据集训练集部分输入到所构建的词嵌入模型中进行特征提取，将测试集用词嵌入模型输出的特征向量矩阵表示，将其通过深度LSTM自编码器联合聚类模型进行进一步聚类及分析，得到较优的聚类性能指标，如图8及表2所示。

表2不同方法下SRR492190测试集聚类结果比较

将通过词嵌入特征向量矩阵表示后的数据直接进行K-means聚类，精确率高达83.0％，召回率为76.9％，调整兰德系数(Adjusted Rand Index，ARI)为80.2％，将数据输入至构建好的深度聚类模型，所得到的最优聚类性能进一步提升，调整兰德系数(AdjustedRand Index，ARI)高达91.3％，精确率高达94.0％，召回率为89.4％。将测试数据集直接通过其他几种对比实验进行聚类，聚类性能与本实施例相比结果较差，各个性能指标均未达到80％，其原因可能为：测试集所选取的数据集较少，其他方法无法通过较少的数据集得到较为准确的特征向量，从而无法实现更好的聚类；本实施例相较其他方法而言，先选取数据集中训练集部分进行特征提取，将特征向量矩阵训练到一定程度，用其表示该数据集的测试集会得到更好的聚类性能，并且本实施例所选取的数据集为不平衡数据集，利用词嵌入特征提取模型与深度聚类模型相结合，进一步改善损失误差函数以及聚类层聚类中心的更新方式，可得到更优的聚类结果。将未知物种通过训练好的特征向量矩阵表示，并通过本实施例构建的深度聚类模型进行聚类，按照结果进行拼接，并计算拼接后得到的重叠群(contigs)完整度及污染度。针对SRR492190数据集而言，未知标签序列共1020842条，read1和read2各包含510421条序列，由于相对应的read1和read2一定属于同一物种，故只需对read1进行深度聚类即可，read2序列的物种分布与read1深度聚类结果一致。由于所选取的宏基因组中已知标签部分为极不平衡数据集，其中最高丰度物种序列占总体的81.65％，故未知物种序列分布与其应大致相似，亦为不平衡数据集。经过深度聚类后，如图9、图10及表3所示：

表3不同方法SRR492190未知物种片段聚类拼接后序列完整度与污染度比较

本实施例共发现了1个未知物种(完整度＞50％，污染度＜6％)，完整度高达74.23％，污染度为4.88％。将未知物种序列通过其他几种方法进行对比实验，均发现了1个物种，其中Sparse Coding(2020)方法组装的重叠群(contigs)污染度为4.40％，但其完整度较低，仅为67.68％；相较于其他几种方法的对比实验而言，利用通过本实施例构建的模型进行深度聚类，按照聚类结果进行拼接的重叠群(contigs)在总体性能指标上具有一定的优势，完整度较高，污染度也相对较低。

综上所述，本发明提供了一种基于预训练和深度聚类的宏基因组物种重建方法。基于宏基因组DNA片段的深度LSTM自编码器联合聚类方法，设计了基于图卷积神经网络联合Focal Loss损失函数的词嵌入特征提取模型以及基于LSTM自编码器联合改进的FCM算法的深度聚类模型。本发明将所选取的宏基因组数据集中每条DNA片段用滑窗法以k-mer的形式编码，应用两层GCN网络作为词嵌入模型训练高维特征，将宏基因组中没有比对标签的片段通过上述模型训练后得到的词嵌入结果表示，作为后续深度聚类的数据样本。本发明构建了一种基于深度LSTM自编码器联合聚类算法模型，对上述处理后的数据集进行聚类分箱，将深度学习与聚类结合在一起，重构误差与聚类误差同步优化，相比于其他算法，可以进一步提升二者性能，同时计算量也较小。最后，利用拼接软件对聚类后的片段按照类别进行组装，并用check-m工具计算组装后的重叠群(contigs)的完整度与污染度。在用户使用时，只需要针对所选取的数据集的大小及序列长度对整个模型的参数进行调整，重新运行模型即可得到聚类结果，大大提高了准确度与便利性。经过真实数据集的验证以及与现有的其他多种方法的聚类结果进行比较，结果证实，本发明能够得到更加优秀的聚类结果。不仅如此，在发现未知物种的层面上，本发明所发现的未知物种相较于其他方法而言，完整度更高，污染度更低。

应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims

1.一种基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，包括步骤：

第一步，特征提取步骤，包括：

针对预处理后的数据集构建词嵌入模型；

构建模型误差函数；

保存输出的特征向量矩阵；

第二步，深度聚类步骤，包括：

计算并分析模型的聚类性能指标；

调整模型参数，得到模型最优的聚类性能；

2.根据权利要求1所述的基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，所述对原始数据集进行预处理具体包括步骤：

a)下载微生物群落的宏基因组序列数据集；

3.根据权利要求2所述的基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，所述低质量序列包括序列中连续存在超过5个碱基N的序列。

4.根据权利要求1所述的基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，所述构建词嵌入模型具体包括步骤：

b)将切割后的序列数据库构建成拓扑图；

c)通过k-mer序列的共现信息来建立两个k-mer节点之间的边；

5.根据权利要求1所述的基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，所述步骤构建模型误差函数中，采用Focal Loss作为特征提取模型的误差函数。

6.根据权利要求1所述的基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，所述编码器和解码器分别包括两层LSTM和两层全连接层，编码器和解码器的中间构建了一层聚类层。

7.根据权利要求6所述的所述的基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，所述LSTM的神经元由若干个递归连接的记忆区块构成，每个记忆区块包括1个记忆单元以及3个逻辑单元，所述逻辑单元包括遗忘门、输入门和输出门。

8.根据权利要求1所述的基于预训练和深度聚类的宏基因组物种重建方法，其特征在于，所述步骤计算并分析模型的聚类性能指标中，所述聚类性能指标包括精确率、召回率以及调整兰德系数。

9.一种存储介质，其特征在于，所述存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如权利要求1-8任意所述的基于预训练和深度聚类的宏基因组物种重建方法中的步骤。

10.一种终端设备，其特征在于，包括处理器，适于实现各指令；以及存储介质，适于存储多条指令，所述指令适于由处理器加载并执行如权利要求1-8任意所述的基于预训练和深度聚类的宏基因组物种重建方法中的步骤。