代谢组数据分析结题报告

温馨提示:请使用火狐或者Chrome的网页浏览器来查看报告

微科盟代谢组数据分析结题报告


一 概述

代谢物亦称中间代谢物,是指通过代谢过程产生或消耗的物质,生物大分子不包括在内。代谢组是指生物体内源性代谢物质的动态整体(代谢物的集合)。而传统的代谢概念既包括生物合成,也包括生物分解,因此理论上代谢物应包括核酸、蛋白质、脂类生物大分子以及其他小分子代谢物质。代谢组目前只涉及相对分子质量约小于1000的小分子代谢物质。代谢组学是对某一生物或细胞在一特定生理时期内所有低分子量代谢产物同时进行定性和定量分析的一门新学科。以组群指标分析为基础,以高通量检测和数据处理为手段,以信息建模与系统整合为目标的系统生物学的一个分支。代谢组研究有非靶向和靶向两种,非靶向代谢组学无偏向地检测样本中所有能检测到的代谢物分子,通过生信方法进行差异分析和通路分析,寻找生物标志物。靶向代谢则是针对特定的代谢物进行检测。

对代谢组进行研究,能够考察生物体系在特定时期受到刺激或扰动前后所有小分子代谢物的组成及其含量变化,从而表征生物体系的整体代谢特征。如:糖、有机酸、脂质、维生素、氨基酸、芳香烃等。也可以 研究代谢物表达量的变化,代谢物与生理病理变化的关系,寻找生物标志物,发现代谢途径。


二 项目流程

2.1 实验流程

详细样品检测步骤请参见实验检测报告。

2.2 数据分析流程

本代谢组分析流程基于R语言MetaboAnalystR包[1]。如图2-2,代谢组的数据分析流程如下:

  1. 样本数据的QC质量控制以及批次校正
  2. 样本数据的标准化
  3. 代谢物含量统计(柱形图,热图)
  4. 无监督的降维分析(PCA)
  5. 特征代谢物(biomarker)的筛选(PLSDA,OPLSDA,单变量分析,火山图,机器学习)
  6. 相关分析
  7. 如果是非靶向代谢组,还可以进行特征代谢物的通路分析(富集分析,拓扑分析,通路图)
图2-2 代谢组数据分析流程


三 数据预处理 (00-DataProcess目录)

3.1 质量评估(QA)和质量控制(QC)(00-DataProcess/01-QA_QC

利用无监督的主成分分析( Principal Component Analysis , PCA ),我们对每一个组的样品建模,然后展示得分图(score plot, 图3-1 )。这种组内的 PCA 分析摒弃了组间的干扰,能让我们更清晰地观察组内的变异并找出可能的离群点( outlier )。图中每个点代表一个样本,而点在图中的位置由该样本中的所有代谢物共同决定。图中的椭圆是基于 Hotelling T2 计算并画出的 95% 置信区间。样本落在椭圆之外暗示该点可能是个离群点。但通常只有在严重离群如将整个模型拉扯得变形时才会考虑删掉该样本。

图3-1 各分组PCA得分图

图中每个点代表一个样本,而点在图中的位置由该样本中的所有代谢物共同决定。图中的椭圆是基于 Hotelling T2 计算并画出的 95% 置信区间。样本落在椭圆之外暗示该点可能是个离群点。

也可利用所有样本点PC1值的分布来评估实验室样本制备、样本测量过程是否处于可控的状态。样本点超出控制界限(3 倍标准差)被认为是离群点。如图3-2, 通常来说,所有点会在控制界限之内。其中,大部分点会在 2 倍标准差之内围绕 平均值 轴上下波动,少量的点会接近控制界限。在大队列的代谢组学研究中,由于仪器的波动,批次间和日间的变异是很常见的。某一个批次的质控样本和测试样本可能会和其他批次的样本有着明显的偏移。因此,在进行后续的统计分析之前,我们可能需要用指控样本来矫正批次间的变异。但是超出三倍标准差的点通常被认为是离群点。如果是不恰当的样本前处理或是进样引起的,离群点样本需要被重新处理和分析。如若不然,这些样本需要进一步的研究,因为可能是真实的生物学差异引起的这些波动。

图3-2 各样本点PC1分布

横坐标为检测样本,纵坐标对应样本在PCA分析中的PC1值。

3.2 标准化校正(Normalization)(00-DataProcess/02-Normlization

在生科云的在线分析中,我们会将是否进行标准化校正的选择权交给用户,然而,我们还是要在这里强调标准化校正的效果和重要性。标准化校正分三步:

  1. 样本内校正:即样本里所有features的丰度除以该样本的丰度中位数(类似相对丰度计算);这是为了校正library size,在测定过程中,各个样本得到的总代谢物含量差异巨大,而这种差异通常是样本采集和测定过程带来的系统误差,除以中位数,平均值,或者总和,都是校正这种系统误差的常用方式;
  2. 含量矩阵校正:即对所有含量值进行log转化;T检验以及ANOVA这样的差异比较方法,都要求代谢物含量服从正态分布,因此,我们一般用log转化,使得代谢物含量分布接近正态分布
  3. feature内校正:即feature对应的所有样本丰度减去该feature丰度均值再除以该feature丰度标准差;feature内校正的目的是使得所有代谢物的均值和标准差(或者说中位数,四分位点,scale)在同一个水平上;PCA,PLSDA,OPLSDA以及机器学习这样的分析,如果不进行代谢物的标准化,那么均值和标准差高的代谢物的重要性将会倾向于高于均值和标准差低的代谢物的重要性,而这样的结果显然不是我们想要的结果,在分组间差异大的才应该是重要性高的。

如图3-3所示,在标准化校正前,代谢物含量的中位数和上下四分位点参差不齐,差异较大,但标准化校正后,基本都在同一个水平上了。

图3-3 标准化校正前后各个代谢物的含量分布

含量分布用箱线图表示,箱线图从左到右分别对应离群点,最小值,下四分位点,中位数,上四分位点,最大值,离群点。左边图为标准化校正前的分布,右边图为标准化校正后的分布。


四 代谢物含量统计(01-ConcentrationSummary目录)

4.1 百分比含量柱形图(01-ConcentrationSummary/1-Barplot)

计算每个样本内各个代谢物的含量百分比,然后用堆积柱形图可视化,能直观的比较分组之间的代谢物组成结构差异。图4-1展示了含量排在前20的代谢物,其余的代谢物被包括进Others中。

图4-1 含量前20的代谢物百分比堆积柱形图

横坐标是样品名,根据分组顺序排序,同时用不同颜色标注了不同的分组样本。纵坐标表示各个代谢物的百分比含量,从上而下对应代谢物的柱的顺序与图例一致。

结果文件夹中也提供了分组均值的百分比含量堆积柱形图,请老师自行查看结果文件夹。

4.2 热图聚类图(01-ConcentrationSummary/2-Heatmap)

为了研究不同样品间的相似性,还可以通过对样品进行聚类分析从而构建样品的聚类。选取关注的代谢物(默认选取代谢物含量排名前30),根据样品的代谢物组成,实现样品聚类,以此考察不同样品或者分组间的相似或差异性;也根据代谢物含量在各样本的分布进行聚类,寻找代谢物或样本的聚集规律。如图4-2,如果不同分组的样本聚类到不同位置,说明分组间的这30种代谢物组成结构差异较大。

图4-2 代谢物热图聚类结果

纵轴为样品名称信息,同时也包括了分组信息。横轴为代谢物。图中上方的聚类树为代谢物在各样本中分布的相似度聚类,左侧的聚类树为样品聚类树,中间的热图是代谢物含量热图,颜色与代谢物含量(Z-Score)的关系见图右上方的刻度尺。

结果文件夹中也有不聚类根据样品分组顺序排序的热图,请老师自行查看结果文件夹。


五 无监督的PCA分析(02-PCA目录)

PCA分析是无监督的,也就是说PCA分析并不考虑样品的分组信息。PCA旨在挑选一些相互独立的能代表所有代谢物的因子(可以理解为所有代谢物的加权和),这些因子要能最大程度地保留原代谢物的变异信息(方差),挑选保留信息最多的两个因子,画散点图,能看出各个样本间的代谢物组成结构差异。如图5-1,如果图中两个分组的点云(point cloud)明显分布在不同区域,说明两个分组的代谢物组成结构差异较大。

图5-1 PCA图

PCA图中,每个点对应一个样本,两个点之间的距离近似与两个样本的代谢物组成结构差异(欧式距离)。不同的分组用不同的颜色标出,椭圆标示的区域是样本点95%置信区域。


六 有监督的最小二乘判别分析

6.1 偏最小二乘判别分析(03-PLSDA目录

偏最小二乘判别分析是有监督的,即分析时需要提供分组信息,PLS-DA寻找能够寻找最大程度区分样品分组的因子(因子可以理解为所有代谢物的加权和)。判别分析将要预测的不连续分类变量编码为一个潜变量,潜变量是连续的,这样就可以在解释变量和潜变量之间建立回归,用最小二乘回归的理论来求解。偏最小二乘判别分析是线性判别分析的升级版,适用于解释变量存在大量共线性问题的组学数据。PLS-DA通过投影分别将预测变量和观测变量投影到一个新空间(新空间各个维度间是相互独立,没有共线性问题的),来寻找一个线性回归模型。如图6-1,我们用区分效果最好的两个因子来画散点图,如果不同分组的样本的点云分布在不同区域,说明PLS-DA模型的判别效果较好,分组间应该存在有显著差异的代谢物。

图6-1 PLS-DA点云图

图中,每个点对应一个样本,横纵坐标是判别效果最好的两个因子的值。不同的分组用不同的颜色标出,椭圆标示的区域是样本点95%置信区域。

有时只看图6-1来判断PLS-DA模型的区分效果(预测能力),不具有说服力。这时候我们可以运用置换检验(Permutation Test)的方法,先构造置换检验统计量(这个统计量通常能衡量模型的预测能力),用置换的方法来求检验统计量的分布(这个分布是随机分布,由随机因素形成的),把样本实际观测到的检验统计量与检验统计量分布进行比较,判断我们的PLS-DA判别模型的判别效果是否是完全随机因素导致的。如图6-2,PLSDA的置换检验选择预测准确度(prediction accuracy)作为检验统计量,如果观测检验统计量在随机分布右侧(观测值明显大于随机值),p值小于0.05,说明我们的PLS-DA判别模型显著,判别效果显著,可以显著区分不同的分组,也即分组间应存在显著差异的代谢物。

图6-2 PLS-DA置换检验的检验统计量分布以及p值

横坐标表示置换检验检验统计量(模型预测准确度)的分布区间,纵坐标是置换过程中该区间的检验统计量的出现频次,箭头所指位置是实际观测到的检验统计量值,如果这个值离随机分布较远,说明模型区分效果不是随机的,模型区分效果显著。

同时也能够寻找在判别分析过程中作用(重要性,VIP, Value Importance in Projection)最大的代谢物,这些代谢物可以作为区分不同分组的biomaker。一般来说,VIP大于1的代谢物对判别分析的贡献是较大的,这些代谢物在分组间有较大差异。如图6-3,黄色区域中标出代谢物名称的代谢物,是校正后p值小于0.05,VIP大于1的代谢物,这些代谢物在分组间有显著差异,在PLSDA中起到重要作用,应该重点关注。

图6-3 PLS-DA代谢物重要性图

图中,每个点代表一个代谢物,横坐标是VIP(Value Importance in Projection),纵坐标是FDR校正后的p值(log10转化)。

6.2 正交偏最小二乘判别分析(04-OPLSDA目录

PLSDA利用代谢组数据所有成分进行预测,往往会造成较为严重的过度拟合(模型在训练样本群体中表现很好,但是在验证样本群体中却表现一般)。正交偏最小二乘判别分析是PLSDA的改进。在正交偏最小二乘判别分析中,回归模型是建立在包含分组信息的解释变量(代谢组数据)和响应变量(分组信息)之间的,模型会将与分组无关的信息过滤掉。与PLSDA相比,OPLSDA的优势在于只利用代谢组数据的一个成分(因子,componet)去预测分组,而其他成分都用来描述与这个成分(预测成分)正交的(无关的)变异[2]。在OPLSDA中,我们通常用R2Y和Q2来衡量模型的区分效果,两个值通常比较接近,Q2=1-模型误差方差/模型总方差,R2Y=模型能解释分组差异的方差/模型总方差,R2Y和Q2越接近1,说明模型的区分效果越好。结果文件夹中表格”oplsda_model_fitness_summary.csv“还提供了R2X,R2X=模型能解释代谢组数据变异的方差/模型总方差,R2X不常用。

OPLSDA点云图和OPLSDA代谢物重要性图与PLS-DA相似,请参照PLS-DA部分进行解读

在OPLSDA的置换检验中,我们用Q2作为检验统计量,用置换的方法求得Q2的随机分布。如图6-4,如果箭头所指的实际观测Q2在随机分布右侧(观测值明显大于随机值),说明Q2不是随机的,是显著的,模型的预测能力显著,也即分组间应该存在显著差异的代谢物。

图6-4 OPLS-DA置换检验的检验统计量(Q2)分布以及p值

分布图是Q2的置换随机分布,箭头所指的是实际观测的模型Q2。


七 单变量分析(05-UnivariateAnalysis目录

某些情况下,我们不仅想知道代谢物在分组间是否发生了变化(差异是否显著),同时也想知道这个变化有多大,以此来评估代谢物变化会给生物体造成多大的影响。我们可以通过计算代谢物变化的倍数(Fold Change, FC)来衡量这个变化的大小,结合p值,来筛选一些重点关注的代谢物。默认情况下,我们会将分组名按字母数字顺序排在前面的分组作为参照,计算另一个分组的均值和参照组的均值的变化倍数,上调变化倍数为正,下调倍数为负。如图7-1,黄色区域的代谢物是p小于0.05,同时变化倍数绝对值大于2的代谢物,这些代谢物在分组间差异显著,且变化较大,应该重点关注。

图7-1 倍数变化火山图

每个点代表一个代谢物,横坐标为变化倍数,纵坐标为T检验p值,变化倍数越大,p值越小(log10(p)越高),点越大。

为了直观地展示代谢物在组之间的差异,我们对单维统计分析获得的排名靠前(P值较小的前25)的代表性的差异代谢物的作箱式图,如下图7-2所示。

图7-2 代谢物差异箱式图

*,**,***分别对应p<0.05, p<0.01, p<0.001


八 机器学习

机器学习也属于判别分析,它是判别分析的延伸。机器学习是一门多领域交叉学科,涉及概率论、统计学、逼近论、凸分析、算法复杂度理论等多门学科。常见的机器学习方法有支持向量机,随机森林,神经网络等。

8.1 随机森林(06-RandomForest目录

随机森林算法是套袋法(bagging)和决策树的结合,大致过程如下:

结果文件夹中,图“RF_outliersdpi300.pdf”中标出了在随机森林分析中,容易判别出错的样本名称,这些样本和组内其他样本差异较大。

在随机森林分析中,我们用“Mean Decrease Accuracy”和“Mean Decrease Gini”来衡量一个代谢物在随机森林中判别分组的重要性。把一个代谢物的取值变为随机数,随机森林预测准确性的降低程度就是“Mean Decrease Accuracy”;“Mean Decrease Gini”则是一个代谢物对分类树所有节点上观测值的异质性的影响。两个值越大,代谢物在随机森林中的重要性越大。图8-1显示了在随机森林中最重要的15种代谢物,这些代谢物应该在分组之间有较大的差异。

图8-1 随机森林中最重要的15种代谢物

左边图的横坐标是“Mean Decrease Accuracy”,它衡量了一个代谢物在随机森林中的重要性;右边图是15中代谢物在两个分组里的含量的热图。

8.2 支持向量机(07-SupportVectorMachine目录

支持向量机(support vector machines, SVM)是一种二分类模型,它的基本模型是定义在特征空间上的间隔最大的线性分类器,间隔最大使它有别于感知机;SVM还包括核技巧,这使它成为实质上的非线性分类器。SVM的的学习策略就是间隔最大化,可形式化为一个求解凸二次规划的问题,也等价于正则化的合页损失函数的最小化问题。SVM的的学习算法就是求解凸二次规划的最优化算法。

SVM模型的最终预测效果可以用ROC曲线来表示,如图8-2。AUC曲线下面积越接近1,说明模型预测效果越好,也说明分组间代谢物的差异越大。

图8-2 SVM ROC 曲线

横坐标为准确度,纵坐标为敏感度,准确度和敏感度越高,模型的预测效果越好。AUC为曲线下面积,曲线下面积越接近1,模型预测效果越好。

同时,我们用套袋法和SVM结合,即用Bootstraping的到k个训练集,在每个训练集里独立挑选最重要的代谢物训练SVM。计算每个代谢物在k个训练集中被选中的频率,这个频率可以用来衡量代谢物在SVM分析中的重要性,同时也可以计算代谢物在所有训练集中的平均重要性。图8-3展示了在SVM分析中最重要的15个代谢物,这些代谢物在分组间应该有较大差异。

图8-3 SVM中最重要的15种代谢物

左边图的横坐标是SVM平均重要性,右边图是15中代谢物在两个分组里的含量的热图。


九 相关分析(08-CorrelationAnalysis目录

通过对比每一个组的代谢物相关性分析谱图,我们可以观察代谢物的相关性在不同组间是否有变化。因为在病理条件下,某些代谢物之间的相关性可能会改变,如从正相关变为不相关或负相关。这将会提示我们关注在这些代谢物所在的通路。我们提供了单个分组以及所有分组一起的pearson相关分析结果的表格,聚类和不聚类的热图,如图9-1是所有分组一起的聚类的相关性热图。

图9-1 Pearson相关性热图

总含量前100的代谢物的相关系数矩阵,相关系数由颜色表示,红色正相关,绿色负相关。


十 参考文献

[1] Jasmine C , Jianguo X , Oliver S . MetaboAnalystR: an R package for flexible and reproducible analysis of metabolomics data[J]. Bioinformatics, 2018.
[2] Westerhuis J A , Velzen E J J V , Hoefsloot H C J , et al. Multivariate paired data analysis: multilevel PLSDA versus OPLSDA[J]. Metabolomics, 2010, 6(1):119-128.


十一 交付结果目录结构

交付数据分为若干个子目录。

|ResultsMetabonome/  [主要的结果文件]
|-- 00-DataProcess / [数据预处理相关信息]
|-- 01-ConcentrationSummary / [代谢物含量统计]
|-- 02-PCA / [PCA分析结果]
|-- 03-PLSDA / [PLSDA分析结果]
|-- 04-OPLSDA / [OPLSDA分析结果]
|-- 05-UnivariateAnalysis / [单变量分析]
|-- 06-RandomForest / [随机森林分析结果]
|-- 07-SupportVectorMachine / [支持向量机]
|-- 08-CorrelationAnalysis / [相关分析结果]
|-- 09-PathwayTopoEnrichment / [通路分析结果,富集分析,拓扑分析](非靶向限定)
|-- FiguresTablesForReport / [本报告的图表文件夹]
|-- 结题报告.html / [结题报告]


十二 联系我们

地址:广东省深圳市南山区南海大道3688号
邮编:518060