微科盟宏基因组和理化组关联分析结题报告

---

一、概述

宏基因组分析内容包括物种组成分析、丰度分析、功能注释等，旨在获取宏基因组的关键特征，例如微生物种类、丰度信息和潜在的功能通路； 理化组学分析内容包括化合物组成分析、浓度分析、差异分析等，旨在获取理化组学的关键特征，例如化合物种类、浓度信息。

宏基因组和理化组关联分析是一种综合性的研究方法，运用多种统计学方法（如相关性分析，偏最小二乘分析等）来揭示宏基因组和理化组之间的关联性，进而确定哪些微生物或功能与理化指标之间显著相关。流程旨在深入理解生物体内微生物群落（宏基因组）和理化指标之间的相互关系，以揭示它们在生物体健康和疾病中的功能和相互作用。

---

二、项目流程

关联分析流程基于R语言多个R包^[1-8]。流程图1，关联分析的具体步骤如下：

1.理化指标数据标准化（样本内校正，含量矩阵校正，feature内校正），宏基因组数据标准化（hellinger标准化）。

2.将关联分析分为两类：理化指标-物种关联分析，理化指标-功能关联分析。

3.相关性分析。

4.回归分析。

5.非监督类模型关联分析：PCA分析, Procrustes分析, PLS分析, sPLS分析, PLS分析, RDA/CCA分析。

6.监督类模型关联分析：随机森林，支持向量机，模型的ROC评估曲线。

---

三、数据预处理

3.1 标准化校正

01.Normalization/
├── Physicochemical          [理化组数据]
|   ├──input_normalization.tsv                     [经过标准化的理化组丰度表]
|   └──input_sample_normalization.tsv              [只经过样本内校正的理化组丰度表]
└── metagenomic              [宏基因组数据]
├── All.KO.abundance_unstratified.hellinger_norm.tsv     [hellinger标准化后的KO丰度表]
└── All.Taxa.OTU.hellinger_norm.tsv                      [hellinger标准化后的物种注释丰度表]

在生科云的关联分析中，我们将理化组数据进行了标准化校正，校正的步骤分为三步：

1.样本内校正：即样本里所有features的丰度除以该样本的丰度中位数（类似相对丰度计算）；这是为了校正library size，在测定过程中，各个样本得到的总理化指标数值差异巨大， 而这种差异通常是样本采集和测定过程带来的系统误差，除以中位数,平均值,或者总和，都是校正这种系统误差的常用方式；

2.含量矩阵校正：即对所有含量值进行log转化；T检验以及ANOVA这样的差异比较方法，都要求理化指标数值服从正态分布，因此，我们一般用log转化，使得数据分布接近正态分布；

3.feature内校正：即feature对应的所有样本丰度减去该feature丰度均值再除以该feature丰度标准差；feature内校正的目的是使得所有理化指标的均值和标准差（或者说中位数，四分位点，scale）在同一个水平上； PCA，PLSDA，OPLSDA以及机器学习这样的分析，如果不进行数据的标准化，那么均值和标准差高的feautres的重要性将会倾向于高于均值和标准差低的features的重要性，而这样的结果显然不是我们想要的结果，在分组间差异大的才应该是重要性高的。

除了理化组数据的标准化步骤，我们使用Hellinger转换方式对宏基因组数据也进行了标准化预处理。标准化处理可以使不同组学数据的尺度一致。理化组数据和宏基因组数据可能具有不同的测量单位和数据范围，通过标准化，可以将它们调整到相同的尺度，使得它们可比性更强。Hellinger转换是宏基因组数据常用的转换方式，特别适用于处理丰度数据。这种转换方式通过将丰度数据的平方根与极端值缩小，降低了高丰度微生物的权重，减轻了不同维度之间的差异对分析结果的影响。

---

四、关联分析内容

4.1单变量关联分析

02.Univariate_Correlation_Analysis/
├──Function.VS.Physicochemical             [功能注释-理化指标关联分析]
|  ├──KEGG
|  |  ├──heatmap  [相关性热图]
|  |  └──network  [相关性网络图]
|  ├──GO
|  ├──EggNOG
|  └──...
└──Species.VS.Physicochemical              [物种注释-理化指标关联分析]
   ├──heatmap  [相关性热图]
   └──network  [相关性网络图]

在理化组和宏基因组的关联分析中，网络图和相关性热图的绘制为研究提供了直观而全面的视觉展示，有助于深入理解两者之间的复杂关系。

4.1.1相关性热图

相关性热图先筛选丰度较高的features，然后应用Spearman相关性分析法进行分析，展示了理化指标和微生物的相关性系数，颜色的深浅表示相关性的强弱；P值小于0.001则标记为***，0.001<P<0.01则标记为**，0.01<P<0.05则标记为*。通过观察整个相关性热图，可以了解理化指标和微生物整体上的关联趋势。深色区域表示强烈的相关性，而浅色区域表示较弱的相关性；进而可以识别出理化指标与微生物之间的模式，例如聚类结构、明显的正相关或负相关群体。这有助于理解它们之间的协同变化和可能的生物学意义。 图中红色表示正相关，蓝色表示负相关；纵轴为宏基因组features(如果是物种注释，则包含物种分类信息)或者功能注释，横轴为理化指标。

4.1.2相关性网络图

网络图展示了理化指标和微生物群落之间的关联关系，每个节点代表一个理化指标或微生物，而连线表示它们之间的关联；网络图能够发现理化指标和微生物之间的正负相关关系，有助于生物学上的解释。通过节点的颜色、形状和大小等信息，可以推测出理化指标与微生物之间可能的共生关系、代谢途径的调控和互补作用。

4.2多变量关联分析

4.2.1非监督类模型分析

4.2.1.1Procrustes分析

03.Multivariate_Correlation_Analysis/
├──Function.VS.Physicochemical             [功能注释-理化指标关联分析]
|  ├──KEGG
|  |  ├──01.Procrustes           [普鲁克分析]
|  |  |  ├──PCA_analysisReport_Omic1.txt       [理化指标PCA分析报告]
|  |  |  ├──PCA_analysisReport_Omic2.txt       [宏基因组PCA分析报告]
|  |  |  ├──PCA_Biplot2.png                    [两个组学的PCA双标图]
|  |  |  ├──Procrustes_analysis_sites.tsv      [普鲁克分析的坐标信息]
|  |  |  ├──Procrustes_residuals.tsv           [普鲁克分析的残差值信息]
|  |  |  ├──procrustes_resultPlot.png          [普鲁克分析图]
|  |  |  ├──Procrustes_resultReport.txt        [普鲁克分析结果报告]
|  |  |  ├──protest_resultReport.txt           [验证结果报告]
|  |  |  └──residuals_plot.png                 [残差值图]
|  |  ├──02.PLS                  [PLS分析]
|  |  ├──03.RDA_CCA              [RDA/CCA分析]
|  |  ├──04.RandomForest         [随机森林分析]
|  |  └──05.SVM                  [支持向量机分析]
|  ├──GO
|  ├──EggNOG
|  └──...
└──Species.VS.Physicochemical              [物种注释-理化指标关联分析]
   └──...

Procrustes分析是一种用于比较两个或多个数据集之间的相似性的方法^[3]， 由于宏基因组与理化组中，存在大量低丰度features，而这些features对于相关性分析贡献度极低，所以进行该模块分析之前，均根据丰度进行了features筛选，Procrustes分析适用于组学关联分析。这种分析通常涉及到PCA双标图、普氏分析图和残差值图等步骤：

PCA双标图：如图4-1通过主成分分析（PCA），将多维数据降维到较低维度， 然后在双标图上展示样本和变量在这些主成分上的分布，以便观察样本之间和变量之间的关系。 其中样本和变量在主成分上的投影用箭头和点表示。样本之间的距离和角度反映它们在主成分方向上的相似性或差异性。 不同箭头的方向之间的角度可以反映原始特征之间的相关性。如果两个箭头的方向相近，表示对应的特征在数据中可能存在相关性。 箭头的长度表示数据在主成分方向上的方差大小。较长的箭头表示在相关主成分上有更大的变化。因此， 箭头的长度可以用来衡量特征的重要性，对应的特征在数据变异性中的贡献更大。

普鲁克分析图：通过普鲁克分析（Procrustes Analysis）对两个数据集进行旋转变换， 使它们在最佳方式下重叠，进一步揭示它们的相似性。结果中得到的旋转矩阵， 将一个数据集映射到另一个数据集的坐标系中。旋转图展示了这种映射关系， 帮助观察样本或变量之间的相似性。每个线段代表一个样本，线段一端实心圆点代表微生物组数据样本点， 线段另一端实心三角形代表相同样本的理化组数据样本点；连线代表两排序构型的矢量残差， 可评价二者间的变异情况，连线越短，表示两个数据集之间一致性越高。M平方代表偏差平方和，是数据离散程度的指标。

残差值图：进行完普鲁克分析之后，通过比较两个数据集在主成分方向上的残差，揭示它们之间的差异。 残差是指样本点到主成分方向上的投影点之间的距离，残差值图展示了这些残差， 可以帮助识别哪些样本在两个数据集之间差异较大。如果样本中两个组学的组成更为相似， 则旋转图中两个匹配的点距离越近，残差较小，反之两个组学无规律，则残差更大。

4.2.1.2 PLS分析

03.Multivariate_Correlation_Analysis/
├──Function.VS.Physicochemical             [功能注释-理化指标关联分析]
|  ├──KEGG
|  |  ├──01.Procrustes           [普鲁克分析]
|  |  ├──02.PLS                  [PLS分析]
|  |  |  ├──ModelBuilding_parameters.csv                [PLS建模参数预测结果]
|  |  |  ├──PLS_Features_correction_BarPlot.png       [载荷值排序图]
|  |  |  ├──PLS_joint_loadings.png                    [联合载荷图]
|  |  |  ├──PLS_joint_scores.png                      [联合得分图]
|  |  |  ├──PLS_loadings.csv                          [联合载荷值信息表]
|  |  |  └──PLS_scores.csv                            [联合得分值信息表]
|  |  ├──03.RDA_CCA              [RDA/CCA分析]
|  |  ├──04.RandomForest         [随机森林分析]
|  |  └──05.SVM                  [支持向量机分析]
|  ├──GO
|  ├──EggNOG
|  └──...
└──Species.VS.Physicochemical             [物种注释-理化指标关联分析]
   └──...

偏最小二乘法（partial least squares，PLS）PLS分析是一种多变量分析方法，主要用于探索两组变量之间的关系，以揭示变量之间复杂的相互作用和关联。PLS通过寻找最能解释X和Y之间共变异的成分来简化模型，有助于在小样本情况下提高模型的稳定性和预测能力。 PLS分析通过对两个数据组间的整合分析，评估两个数据集之间的内在相关性^[4]。PLS模型一方面可反映不同数据组间的整体影响， 另一方面可直接体现不同变量在模型中的权重，权重越大，意味着该变量的变化对另一个组学的扰动更剧烈。该分析模块分为三个部分： 联合载荷图，联合得分图，关联程度最大的features柱状图。

联合载荷图（Joint loadings plot）：该图显示了两个数据矩阵之间的共同结构，即它们之间的关联性。 如果该图是二维散点图，则X轴和Y轴分别代表了PLS第一成分loadings值和PLS第二成分loadings值，每个点表示一个features，点的位置和方向表示该features在两个数据集中的权重和方向。联合载荷图的X轴和Y轴值的大小和方向反映了变量在两个组学数据集之间的关系； 正值和负值表示features在两个数据集中的相关性方向，值的大小表示features对另一个组学的关联程度。features在图中的位置和包裹趋势则表明它们在两个数据集中的共同结构上的重要性和方向。 这些观察可以帮助解释两个组学数据集之间的关联性和共同模式。

注：(a)，物种注释-理化指标 联合载荷图；(b)，功能注释-理化指标 联合载荷图

联合得分图（Joint score plot）：Joint scores图展示了每个样本在PLS中的投影，即它们在共同结构上的表现。 这反映了样本在两个数据集之间的关联性。正负值反映了样本在两个组学数据集之间的相关性方向； 样本在图中的相对位置反映了它们在两个组学数据集之间的关联性， 相邻或聚集的样本可能在两个数据集中表现相似。

注：(a)，物种注释-理化指标 联合得分图；(b)，功能注释-理化指标 联合得分图

最后选择前两个维度（joint loadings1、joint loadings2）载荷值长度前15的features（关联程度最大）绘制柱状图。

注：(a)，物种注释-理化指标 载荷值排序图；(b)，功能注释-理化指标 载荷值排序图

4.2.1.3RDA/CCA分析

03.Multivariate_Correlation_Analysis/
├──Function.VS.Physicochemical             [功能注释-理化指标关联分析]
|  ├──KEGG
|  |  ├──01.Procrustes           [普鲁克分析]
|  |  ├──02.PLS                  [PLS分析]
|  |  ├──03.RDA_CCA              [RDA/CCA分析]
|  |  |  ├──DCA_output.xls                                      [DCA分析结果表格]
|  |  |  ├──groups_RDA.envfit.xls                               [理化指标对宏基因组影响值表格]
|  |  |  ├──groups_RDA.Factors.PERMANOVA.xls                    [PERMANOVA分析报告_features因子]
|  |  |  ├──groups_RDA.features.xls                             [features影响值表格]
|  |  |  ├──groups_RDA.Group.PERMANOVA.xls                      [PERMANOVA分析报告_样本分组]
|  |  |  ├──groups_RDA.sample.xls                               [样本影响值表格]
|  |  |  ├──groups_RDA_features_location_plot.png               [features RDA双序图]
|  |  |  ├──groups_RDA_sample_location_plot.png                 [样本 RDA双序图]
|  |  |  └──groups_RDA_sample_location_plot_with_labels.png     [样本 RDA双序图(带标签)]
|  |  ├──04.RandomForest         [随机森林分析]
|  |  └──05.SVM                  [支持向量机分析]
|  ├──GO
|  ├──EggNOG
|  └──...
└──Species.VS.Physicochemical               [物种注释-理化指标关联分析]
   └──...

RDA和CCA都属于冗余分析（Redundancy Analysis）的范畴， 都是在解释一个变量组对另一个变量组的方差贡献时进行的， 云流程会根据数据特点自动挑选RDA或者CCA分析。 旨在用统计学手段建立两个组学features之间的回归关系，保留能解释的变异， 后将features进行降维，通过夹角指示相关性。箭头的长度表示理化组在RDA轴上的投影的强度， 较长的箭头表示该理化指标在相应RDA轴方向上的贡献较大，与宏基因组的变化有关， 长箭头指示了理化指标与宏基因组数据的关联性更强，对相应的轴有更大的贡献； 箭头与宏基因组features点的夹角表示理化指标与宏基因组的关联性，以及这个关联性的方向。 小夹角表示理化指标和相应的宏基因组特征在RDA轴上有相似的变化趋势， 表明它们之间可能存在关联性。 样本RDA双序图中的p值反映了分组对微生物差异的解释程度，p小于0.05,则解释方差显著。

注：(a)，物种注释-理化指标 features RDA双序图；(b)，功能注释-理化指标 features RDA双序图

注：(a)，物种注释-理化指标 样本RDA双序图；(b)，功能注释-理化指标 样本RDA双序图

4.2.2监督类模型分析

4.2.2.1随机森林

03.Multivariate_Correlation_Analysis/
├──Function.VS.Physicochemical             [功能注释-理化指标关联分析]
|  ├──KEGG
|  |  ├──01.Procrustes           [普鲁克分析]
|  |  ├──02.PLS                  [PLS分析]
|  |  ├──03.RDA_CCA              [RDA/CCA分析]
|  |  ├──04.RandomForest         [随机森林分析]
|  |  |  ├──Physicochemical_model_Error_Rate_Curve.png                                   [理化组RandomForest模型错误率曲线图]
|  |  |  ├──Physicochemical_model_Probability_Graph_Of_Correctly_Classified.png          [理化组RandomForest模型样本裕度图]
|  |  |  ├──Physicochemical_RFimportance.tsv                                             [理化组RandomForest模型中features重要性信息]
|  |  |  ├──PhysicochemicalROC_curve_specificities_ofRandomForest.tsv                    [理化组RandomForest模型ROC曲线信息]
|  |  |  ├──Metagenomic_model_Error_Rate_Curve.png                                   [宏基因组RandomForest模型错误率曲线图]
|  |  |  ├──Metagenomic_model_Probability_Graph_Of_Correctly_Classified.png          [宏基因组RandomForest模型样本裕度图]
|  |  |  ├──Metagenomic_RFimportance.tsv                                             [宏基因组RandomForest模型中features重要性信息]]
|  |  |  ├──MetagenomicROC_curve_specificities_ofRandomForest.tsv                    [宏基因组RandomForest模型ROC曲线信息]
|  |  |  ├──RFMeanDecreaseAccuracy.png                                               [两个组学RandomForest模型features重要性图]
|  |  |  └──ROC_of_RFmodel.png                                                       [两个组学RandomForest模型ROC评估曲线图]
|  |  └──05.SVM                  [支持向量机分析]
|  ├──GO
|  ├──EggNOG
|  └──...
└──Species.VS.Physicochemical               [物种注释-理化指标关联分析]
   └──...

随机森林是一种集成学习方法^[1]，它通过构建多个决策树，每个树都在不同的子样本和特征上训练， 最后将它们的结果综合起来。这种集成方法有效地减小了过拟合的风险，并且对于高维数据集表现良好。 随机森林能够输出每个特征的重要性评分，即特征对于模型性能的贡献程度。这些重要性分数可以用来筛选最相关的特征， 帮助识别与实验组和对照组差异相关的生物学特征。随机森林关联分析的目的是利用多个决策树的集成效果，获取两种组学间， 各features的重要性比较，并通过不同组学的单独建模的ROC曲线（受试者回归曲线， Receiver Operating Characteristic curve）对比，评估哪种组学能更好地分离对照组和实验组。

注意：在多分类模型中（分组方案中存在多个分组），绘制ROC曲线通常是通过将每个类别与其他所有类别结合起来进行二分类来完成的。这是因为ROC曲线是针对二分类问题设计的，它显示了真正例率（True Positive Rate，也称为灵敏度）和假正例率（False Positive Rate）之间的关系。 在多分类情况下，可以将每个类别视为一个正例，其他所有类别合并为一个负例，然后分别计算每个类别的真正例率和假正例率，最终将这些结果绘制在同一张ROC曲线上。

4.2.2.2支持向量机

03.Multivariate_Correlation_Analysis/
├──Function.VS.Physicochemical             [功能注释-理化指标关联分析]
|  ├──KEGG
|  |  ├──01.Procrustes           [普鲁克分析]
|  |  ├──02.PLS                [PLS分析]
|  |  ├──03.RDA_CCA              [RDA/CCA分析]
|  |  ├──04.RandomForest         [随机森林分析]
|  |  └──05.SVM                  [支持向量机分析]
|  |  |  ├──Physicochemical_SVM_features_importance.tsv               [理化组SVM模型features重要性信息]
|  |  |  ├──PhysicochemicalROC_curve_specificities_ofSVM.tsv          [理化组SVM模型ROC曲线信息]
|  |  |  ├──Metagenomic_SVM_features_importance.tsv               [宏基因组SVM模型features重要性信息]
|  |  |  ├──MetagenomicROC_curve_specificities_ofSVM.tsv          [宏基因组SVM模型ROC曲线信息]
|  |  |  ├──ROC_of_SVMmodel.png                                   [两个组学SVM模型ROC曲线图]
|  |  |  └──SVM_FeaturesImportance.png                            [两个组学SVM模型features重要性图]
|  ├──GO
|  ├──EggNOG
|  └──...
└──Species.VS.Physicochemical               [物种注释-理化指标关联分析]
   └──...

支持向量机是一种监督学习算法^[5]，其目标是找到一个超平面，将不同类别的样本分隔开。 在生物信息学中，SVM可以用于分类和回归任务。SVM在分类问题中对样本之间的复杂关系有较好的拟合能力。 支持向量机的支持向量也可以提供有关差异性特征的信息。 除了两个组学features比较以外，与随机森林一样，提供组学分别建模的ROC曲线比较。

宏基因组和理化指标的数据通常是高维度的，其中包含了大量的features。随机森林在高维数据中表现优异， 因为它可以自动选择子集进行训练，从而减小过拟合风险。相比之下，支持向量机在高维数据中的性能高度依赖于核函数的选择和参数调整； 随机森林相对于异常值具有较好的鲁棒性，因为它基于多个决策树的综合结果，异常值对于整体的影响较小。支持向量机在某些情况下可能对异常值较为敏感， 因为它试图找到一个全局最优的超平面，受到极端值的影响较大。

注意：在多分类模型中（分组方案中存在多个分组），绘制ROC曲线通常是通过将每个类别与其他所有类别结合起来进行二分类来完成的。这是因为ROC曲线是针对二分类问题设计的，它显示了真正例率（True Positive Rate，也称为灵敏度）和假正例率（False Positive Rate）之间的关系。 在多分类情况下，可以将每个类别视为一个正例，其他所有类别合并为一个负例，然后分别计算每个类别的真正例率和假正例率，最终将这些结果绘制在同一张ROC曲线上。

---

五、分析所用软件的版本

软件	版本
ggplot2	4.2.3
vegan	4.2.3
OmicsPLS	4.2.3
igraph	4.2.3
pheatmap	4.2.3
randomForest	4.2.2
e1071	4.2.3
pROC	4.2.3
pathview	4.2.1
DESeq2	1.26.0
edgeR	3.28.1
mixOmics	6.20.0

---

六、参考文献

[1] Breiman L. Random Forests. Machine Learning, 2001 45(1), 5-32.

[2] Ito K., Murphy D. Application of ggplot2 to Pharmacometric Graphics. CPT Pharmacomet. Syst. Pharmacol, 2013 2, 79.

[3] Rohlf J. Shape Statistics: Procrustes Superimpositions and Tangent Spaces. Journal of Classification, 1999 16(2), 197-223.

[4] Trygg J., Wold, S. Orthogonal Projections to Latent Structures (O-PLS). Journal of Chemometrics: A Journal of the Chemometrics Society, 2003 17(3), 119-128.

[5] Cortes C., Vapnik V. Support-Vector Networks. Machine Learning, 1995 20(3), 273-297.

[6] Love I., Huber W., Anders S. Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with Deseq2. Genome Biology, 2014 15(12), 550.

[7] Smyth K. edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics, 2010. 26(1), 139.