基于机器学习的微生物溶解有机碳含量估测

马云鹏 ， 朱静 ， 崔兴华

新疆农业大学计算机与信息工程学院，乌鲁木齐 830052

微生物群落在不同的生态系统中发挥重要的作用，虽然微生物群系控制着许多宏观性质，但阐明特定微生物与生态系统功能之间的关系仍然是生态学中一个复杂的问题。第二代DNA 测序技术降低了微生物群落中宏基因组数据的获取难度，使得利用微生物探究宏观性质变化成为可能。然而，微生物数据的分析也存在一定难度，每个研究只采用了少量样本，而且产生的测序结果通常都是高维、稀疏的数据［1］。

随着人工智能的发展，机器学习技术提供了一种分析高维数据的新方法［2-3］，并可用于阐明微生物菌群与环境属性及宿主表型之间的关系。研究显示，利用机器学习的随机森林模型和使用微生物菌群数据可以有效鉴定栖息地、宿主和疾病状态，并证明机器学习模型分类准确性较高［4］；Zeller 等［5］利用线性回归模型结合微生物组数据，可以有效区分患者；Ning 等［6］利用支持向量机模型可以有效对口腔微生物群落进行分类；Lo 等［7］利用卷积神经网络结合微生物组数据可以预测宿主的健康状态；Bokulich等［8］认为机器学习结合微生物组数据不仅可以用于分类问题，还可以用于回归问题，相关机器学习算法也可用于筛选微生物组数据重要特征，揭示相关微生物组数据与目标变量的重要程度，并从机器学习的角度为了解微生物组产生的影响提供一种新的视角。

目前大部分的微生物组数据研究主要是探究宿主表型预测，以及微生物自身及相互之间产生的影响。最近的研究表明，土壤中含有大量微生物［9］，增温、森林植被、有机肥等均会影响土壤中的微生物数量，进而影响溶解有机碳（dissolved organic carbon，DOC）含量［10-13］。研究证明细菌微生物群落在碳循环中发挥重要的作用，并对土壤中DOC 的丰度产生重要影响［14］，且第二代测序技术应用的普及为利用细菌微生物操作分类单元（operational taxonomic units，OTU）丰度数据预测所产生的溶解有机碳提供了可能。OTU 丰度单元是最常用的依据微生物组数据利用机器学习预测目标情况的输入特征，在生物信息学中为了便于对菌群进行分析，需要人为地对原始16S rRNA基因数据进行聚类操作。根据不同的相似度水平，对所有的测序数据进行OTU 划分，一般情况下序列间相似性＞97%可以认定为同一个种，即一个OTU 都是同属一个种的微生物。但是微生物组数据样本的数量一般远小于OTU 丰度数据，造成数据高维的特点，传统的统计学方法难以有效地对高维数据进行处理，而机器学习可以从复杂、高维、大量的数据中进行学习，进而用于阐明微生物菌群与目标属性之间的联系［15］。

本研究应用常见的12 种机器学习模型，根据细菌微生物菌群丰度预测来自植物凋落物分解的DOC 含量，其中细菌微生物菌群的丰度被视为模型特征/自变量，DOC 含量被视为模型的目标/因变量。所有数据用于机器学习，选择预测精度较优的前3 种机器学习模型作为基模型，利用包装法、嵌套法及包装-嵌套融合方法进行特征选择，选取对DOC 含量影响最大的特征子集，基于筛选出来的最优特征子集构建溶解有机碳估测模型，发现梯度提升决策树（gradient boosting decision tree，GBDT）模型具有最高精度，并对模型参数进行网格参数搜索优化，建立了DOC 估测模型，以期为利用细菌微生物估测DOC含量提供新方法。

1 GBDT建模与方法

1.1 GBDT建模

GBDT 算法是一种使用分类与回归树（calssification and regression tree，CART）作为决策树的Boosting 集成学习算法，既可以处理回归问题也可以处理分类问题，是利用加法模型和前向分布算法实现学习的优化过程，即GBDT 是一系列CART 回归树的加法组合，通过新建一棵树来拟合之前预测值与目标值的“残差”，即Boosting 集成学习的核心思想为将一系列弱学习器进行组合，构成一个强学习器。GBDT建模如下。

式中，N为模型训练数据的样本数；yi˙为模型输出值；c为使得损失函数达到最小的常数。

令迭代次数为m=1，2，…，M

①对样本i=1，2，…，N，计算损失函数在负梯度的值，作为残差估计。

式中，yi˙为模型输出值；f(xi)为迭代次数m次的上一次预测值。

②对rmi拟合一个回归树，得到第m棵树的叶节点区域Rmi，j=1，2，…，J，用以拟合残差近似值。

③对j=1，2，…，J，利用线性搜索估计叶节点区域的值，使损失函数极小化。

式中，yi˙为模型输出值；fm-1(xi)为迭代次数m次的上一次第i个样本的预测值；γm为迭代次数m次时的节点残差拟合值。

④更新回归树，计算公式如公式（4）。

式中，fm-1(x) 为迭代m-1 次的预测值；γmj为迭代m次第j个叶节点的残差拟合值；I为指示函数。

最终输出模型，见公式（5）。

式中，fM(x)为迭代M次最终的预测模型。

1.2 数据与方法

1.2.1数据预处理 本研究使用的数据来源于Johansen 等［16］研究土壤微生物在植物凋落物分解产生的碳流量作用相关实验，数据已存入NCBI序列档案（序列号：SRP151768）。在Johansen 等［16］的研究中，将收集到的土壤样本中的微生物群落，接种到3 个含有无菌沙粒和松木凋落物的相同微环境中，并在25 ℃下孵育44 d。第44 天，测定微环境中DOC 含量，从相同微环境的一个子集中提取DNA，并在Illumina MiSeq 上对16S rRNA 基因扩增子进行测序，通过对308 个样品的细菌群落谱进行细化，得到了1 023 条序列，共2 521 个OTU丰度（表1）。

表1 OTU部分样表Table 1 OTU partial sample table

首先对原始数据进行数据清洗。原始OTU数据具有高维、稀疏的特征，也造就了大量空值，将原始数据中所有样本共有的OTU 丰度单元数据为0 的列进行删除，删除空缺列后得到1 709 个特定分类单元的丰度数据，因在原始的数据中OTU 丰度有的高达上万，有的为零，很容易出现“大数吃小数”的问题，在机器学习中算法模型容易忽略丰度较小的OTU 特征，因此需要将原始数据进行归一化处理，见公式（6）。

式中，x*为经过归一化后得到的数据；x为特征数据；xmin为特征数据的最小值；xmax为特征数据的最大值。将数据进行归一化处理之后可以消除量级的影响，使结果映射到［0，1］之间，进行归一化后的数据用于机器学习训练，可以有效提高模型的收敛速度及预测精度。

1.2.2特征选择 本研究中主要使用了Embedded 嵌入法中的特征重要性度量（feature importance select，FIS）、Wrapper 包装法中的递归消除特征（recursive feature elimination，RFE）及包装法-嵌入法融合这3种特征选择方法。

嵌入法首先使用基学习模型进行训练，利用基学习模型得到特征的权值系数，然后根据权值系数从大到小的选择特征。嵌入法是一种可以自己选择使用哪些特征的方法，即算法训练和筛选特征同时进行（图1）。

图1 嵌入法特征选择过程Fig. 1 Feature selection process of embedding method

包装法与嵌入法类似，也是通过自身算法对特征进行选择，通过目标函数来获得特征的重要性，并在每次迭代时保留最佳特征或者剔除最差特征，并在下一次迭代时使用以前未使用过的特征来构建模型，直至将所有特征用尽，然后根据保留或剔除特征的顺序来选择一个最佳特征子集（图2）。

图2 包装法特征选择过程Fig. 2 Feature selection process of wrapper method

嵌入法考虑单个特征对目标的影响，模型简单，计算开销比较小，包装法考虑特征子集对目标的影响程度，弥补了嵌入法针对特征筛选不足的情况，但是反复构建模型计算开销比较大。因此本研究首先使用嵌入法筛选出对目标变量影响比较大的特征，然后利用包装法对嵌入法筛选出来的特征进行二次特征筛选（图3）。

图3 包装-嵌套融合法特征选择过程Fig. 3 Feature selection process of wrapper-embedded fusion method

图4 10折交叉验证示意图Fig. 4 10-fold cross-validation schematic diagram

1.2.3交叉验证 交叉验证法（cross validation）是将数据集D 均匀的划为k份且互斥的子集，循环每次用k-1 个子集的并集当作训练集，余下的子集作为测试集，即可进行k次训练和测试。交叉验证相比较原始的一次性划分数据有诸多优点，如可以明显降低一次划分所造成的随机性；模型进行k次训练可以充分的利用原始数据并且提高模型的泛化能力。

1.2.4评价指标 评价指标采用均方根误差（root mean squared error，RMSE）、平均绝对误差（mean absolute error，MAE）及线性拟合优度（coefficient of determination，R2）来评价本研究的模型。

RMSE 为预测值和真实值之间误差平方和的平均数的平方根，用来衡量预测值同真值之间的偏差，见公式（7）。

式中，y为真实值；y^ 为预测值；m为数据量；yi为第i条数据的DOC 实测值；y^i为第i条数据的DOC预测值。

MAE表示预测值和真实值之间绝对误差的平均值，能反映预测值误差的实际情况，见公式（8）。

式中，y为真实值；y^ 为预测值；m为数据量；yi为第i条数据的DOC 实测值；y^i为第i条数据的DOC预测值。

R2的取值范围为［0，1］之间，越接近于1，说明模型的拟合效果越好，见公式（9）。

式中，m为数据量；yi为第i条数据的DOC 实测值；y^i为第i条数据的DOC 预测值，yˉ为m条数据的DOC平均值。

1.2.5实验方案 ①将人工剔除空值的OTU 丰度数据（共1 709 个特征）归一化后应用于初始的12 种机器学习方法。将308 个样本按照10 倍交叉验证法分成训练集和测试集，训练集用于机器学习训练，测试集用于评估机器学习模型，由于数据量小而且特征多，因此为了避免过拟合现象的发生，使用10 倍交叉验证的方式应用于原始数据进行模型参数调参，可以检测模型的泛化性和稳定性，让模型的输出结果更为可靠。如果只将数据进行简单的人为划分，则会增加样本划分的偶然性，使得预测模型变得不准确。将数据进行10折划分后，所有的数据都可用于模型训练，增加了数据利用率，同时通过10 折交叉验证来确定模型的参数，是经过10 次比较得出的相对最优模型，可以一定程度上防止过拟合和欠拟合现象的发生［17］。所有机器学习模型构建使用Python 语言，应用Pycharm 软件结合sklearn 库，各机器学习模型使用默认参数。②将初始的12 种机器学习方法，依照RMSE、MAE、R2这3 种评价指标选择前3 种较优的机器学习模型作为基学习模型用于特征选择，结合包装法、嵌入法、包装-嵌入融合方法选取最优的特征子集。③将选取的最优特征子集用于初始训练表现最好的机器学习模型，使用网格搜索（grid search，GS）参数优化方法进行参数寻优，从而确定预测模型的最优参数。

2 结果与分析

2.1 机器学习预测结果

由表2 可知，应用12 种机器学习模型进行预测，GBDT 表现最优，选取较优的3 种模型，即GBDT、随机森林（random forest，RF）和极限树（extratree，ET）作为基模型用于包装法的递归消除特征（recursive feature elimination，RFE），选择GBDT用于嵌入法的特征重要性度量（feature importance select，FIS）以及嵌入-包装融合法（RFE-FIS）作为对比。

表2 多种机器学习模型预测结果Table 2 Prediction results of multiple machine learning models

2.2 基模型不同的递归消除特征结果

以GBDT 为基模型的递归消除特征得到112个OTU 丰度单元评分最高（图5）。以RF 为基模型的递归消除特征得到88 个OTU 丰度单元评分最高（图6）。以ET 为基模型的递归消除特征得到141 个OTU 丰度单元评分最高（图7）。以上结果表明基模型不同则通过递归消除特征选取的最佳特征子集不同，应用递归消除特征来进行特征选择时，基模型选择至关重要。

图5 RFE（GBDT）特征选择Fig. 5 RFE （GBDT） feature selection

图6 RFE（RF）特征选择Fig. 6 RFE （RF） feature selection

图7 RFE（ET）特征选择Fig. 7 RFE （ET） feature selection

2.3 融合模型选择特征结果

使用GBDT 作为基模型用于特征重要性选择［FIS（GBDT）］得到246 个OTU 丰度单元特征最为重要，然后使用递归消除特征以GBDT 为基模型进行二次特征选择得到140个OTU丰度单元评分最高（图8）。筛选得到的140 个菌群操作分类单元具体编号如表3所示。结果表明，经过特征重要性度量筛选后，对DOC有影响的特征有246个，但是特征重要性度量只考虑了单个特征对DOC的影响，未考虑多特征对DOC的影响，因此使用递归消除特征进行二次特征选择，可对特征进一步缩减，得到了包含140 个特征的最优特征子集。图8 也表明140个特征子集的模型交叉验证得分最高。

图8 RFE-FIS（GBDT）特征选择Fig. 8 RFE-FIS （GBDT） feature select

表3 REF-FIS（GBDT）特征选择OTU表Table 3 RFE-FIS （GBDT） feature selection OTU table

2.4 预测结果分析

将筛选之后得到的特征用于较优的GBDT、RF 及ET 模型可得如表4 所示的预测结果。进一步分析表明，GBDT 算法模型应用经过二次特征选择的140 个特征取得较优的预测结果，利用网格参数搜索法进一步提升模型的精度，网格搜索法在GBDT 模型的参数空间中对参数值进行分割，化为网格，以尝试每种参数，获得模型的最优参数解。sklearn 库中GBDT 算法模型的默认参数学习率（learning_rate）为0.1，构造回归树的数量（n_estimators）为100，最大数深（max_depth）为3，对模型影响较大的3 种参数使用网格搜索如表5所示。经过网格搜索得到GBDT 算法模型的最优参数为学习率0.07，树的数量574，最大深度3，优化前后的对比精度如表6 所示，优化后预测结果与真实值对比如图9 所示，经历过特征选择和参数优化后的GBDT 模型预测值与真实值的线性拟合效果如图10 所示。以上结果表明，经历过特征选择及参数优化的算法模型对比初始模型的预测精度明显提高。根据测试数据可知，在3 种评价指标下进行过数据特征选择以及模型参数优化后的GBDT 模型达到了最优的预测精度，相比应用最初原始的1 709个特征，经过二次特征选择及参数优化后的GBDT 模型精度得到显著提高，筛选出了对DOC 影响至关重要的140 个OTU 丰度特征，降低了特征数，使得该模型的预测值最为接近实际值。

图9 预测值与真实值对比图Fig. 9 Comparison figure between predicted and true values

图10 线性拟合效果图Fig. 10 Linear fitting effect draw

表4 模型预测结果Table 4 Model prediction results

表5 模型参数网格搜索范围Table 5 Model parameter grid search range

表6 参数优化后精度对比Table 6 Precision comparison after parameter optimization

3 讨论

Johansen 等［16］的研究证明，细菌群落相比真菌群落与溶解有机碳之间有更强的联系；Thompson 等［18］利用Johansen 等［16］的数据通过随机森林的特征重要性度量、人工神经网络以及指示种分析3 种特征选择方法，得到共有的86 个OTU 特征。应用随机森林结合筛选得到的86 个特征预测DOC 的含量，预测值与真实值的皮尔逊相关系数为0.636，线性拟合优度为0.404 5，应用人工神经网络得到的预测值与真实值的皮尔逊相关系数为0.676，线性拟合优度为0.456 9。本研究的结果与Thompson 等［18］的结果相似，证明机器学习可以作为高维稀疏微生物数据的特征选择工具，以及利用细菌微生物菌群反推溶解有机碳含量的一种技术手段，并取得了较好的效果，也从侧面印证了Johansen等［16］研究结果的科学性。

本研究与先前研究不同点在于：通过对原始数据进行数据处理和特征选择，并构建了12 种机器学习模型进行对比分析，超过了Thompson 等［18］研究的2种机器学习方法；使用了5种特征选择方法，多于Thompson 等［18］使用的3 种方法；Thompson 等［18］将数据特征降低到86 个，本研究将数据特征降低到140 个，经过参数优化后的GBDT 模型表现最优，DOC 预测值与真实值的皮尔逊相关系数为0.812，线性拟合优度为0.659 9，优于Thompson 等［18］的研究；Thompson 等［18］的研究选择了257个样本作为训练集，51个样本作为测试集，这样简单的划分具有较大的随意性，本研究为了减少样本划分对实验模型的影响，同时降低模型过拟合的风险使用了十折交叉验证，使得结果更具有说服力。

通过网格搜索对GBDT 模型进行参数优化，最优参数为学习率0.07，树的数量574，数模型的最大深度3，此时模型10 折交叉验证的均方根误差为1.722 0，平均绝对误差为1.293 4，线性回归拟合优度为0.659 9，模型预测精度比参数优化前明显提高。从机器学习的角度解释了细菌微生物菌群中特定的140 个操作分类单元与溶解有机碳有较强的关系，证明仅依靠140 个OTU 丰度数据即可较好地得到溶解有机碳的估测值，该模型在一定程度上解决了数据特征问题，而且也简化了研究人员的实验环节，但是该模型也只是估测，无法达到Johansen 等［16］真实实验所得到的数据，所以与真实结果略有不同，但也为利用细菌微生物估测溶解有机碳提供了一种新的方法。

本研究对细菌微生物菌群与溶解有机碳进行了初步探索，但仍存在许多重要的问题值得进一步研究，如本实验的数据是在实验室人为设定条件下得到的测量值，无法真实模拟野外环境；虽然样本数量足够微生物学家研究使用，但是机器学习使用还略有不足；研究数据仅考虑细菌丰度对溶解有机碳的影响，未考虑时间进程和细菌丰度之间综合影响溶解有机碳含量的关系等。未来研究需进行更多实地研究以更好地了解野外环境下的细菌微生物与溶解有机碳之间的关系，以佐证本实验结果。随着实验样本量的增加，新兴的深度学习可以更好地拟合大数据，发现更深层次的特征，取得比机器学习更好的预测效果［19］，使得未来依据细菌微生物来推测一段时间内溶解有机碳的变化趋势成为一种可能。

综上所述，未来机器学习将会成为分析微生物菌群与其所产生的影响的一个重要工具，研究者应积极利用新技术，探索两者之间的结合点，并借助计算机强大的计算能力，使其成为未来微生物菌群分析的重要工具。