代谢组学数据分析的统计学方法*

柯朝甫张 涛武晓岩李 康Δ

代谢组学数据分析的统计学方法*

柯朝甫1张 涛2武晓岩1李 康1Δ

代谢组学是近年发展快速的一门学科，目前在医学、植物学、微生物学、毒理学、药物研发等诸多领域中得到了广泛的应用［1-5］。代谢组学研究产生大量的数据，这些数据具有高维、小样本、高噪声等复杂特征。如何从复杂的代谢组学数据中提取出有价值的信息，筛选出潜在的生物标志物成为近年来代谢组学研究的热点和难点。据此，本文针对目前代谢组学数据分析中的常用统计学方法及其研究进展进行介绍。

代谢组学数据的特点

代谢组学是系统生物学领域中继基因组学和蛋白质组学之后新近发展起来的一门学科，它通过检测生物体在受到外源刺激或基因修饰后其体内代谢物质的变化来探索整个生物体的代谢机制［6］。其研究对象为生物体内所有内源性小分子代谢物（分子量 ＜1000Da），研究手段为高通量检测技术和数据处理方法，最终目标是数据建模和生物标志物的筛选。生物样品如血浆、尿液、组织等，经过GC/MS、NMR、LC/MS等高通量仪器检测后，得到大量的图谱数据，使用XCMS［7］等软件对这些图谱数据进行转换，获得用于统计分析的标准格式的数据。归纳起来，代谢组学数据具有以下特点：

（1）高噪声：生物体内含有大量维持自身正常功能的内源性小分子，具有特定研究意义的生物标志物只是其中很少一部分，绝大部分代谢物和研究目的无关。

（2）高维、小样本：代谢物的数目远大于样品个数，不适合使用传统的统计学方法进行分析，多变量分析容易出现过拟合和维数灾难问题［8］。

（3）高变异性：一是不同代谢物质的理化性质差异巨大，其浓度含量动态范围宽达7～9个数量级［9］，二是生物个体间存在各种来源的变异，如年龄、性别都可能影响代谢产物的变化，三是仪器测量受各种因素影响，容易出现随机测量误差和系统误差，这使得识别有重要作用的生物标志物可能极其困难。

（4）相互作用关系复杂：各种代谢物质可能不仅具有简单的相加效应，而且可能具有交互作用，从而增加了识别这些具有复杂关系的生物标志物的难度。

（5）相关性和冗余性：各种代谢物并非独立存在，而是相互之间具有不同程度的相关性，同时由于碎片、加合物和同位素的存在使得数据结构存在很大的冗余性，这就需要采用合理的统计分析策略来揭示隐藏其中的复杂数据关系。

（6）分布的不规则和稀疏性：代谢组学数据分布不规则，而且数据具有稀疏性（即有很多值为零），因此，传统的一些线性和参数分析方法此时可能失效。

数据的预处理

代谢组学数据分析的目的是希望从中挖掘出生物相关信息，然而，代谢组学数据的变异来源很多，不仅包括生物变异，还包括环境影响和操作性误差等方面。处理手段主要包括归一化（standardization）、标准化（normalization），即中心化（centering）和尺度化（scaling），以及数据转换（transformation）［10］。归一化是针对样品的操作，由于生物个体间较大的代谢物浓度差异或样品采集过程中的差异（如取不同时间的尿样），为了消除或减轻这种不均一性，一般使用代谢物的相对浓度，即每个代谢物除以样品的总浓度，以此来校正个体差异或其他因素对代谢物绝对浓度的影响。标准化是对不同样品代谢物的操作，即统计学意义上的变量标准化。标准化的目的是消除不同代谢物浓度数量级的差别，但同时也可能会过分夸大低浓度组分的重要性，即低浓度代谢物的变异系数可能更大。数据转换是指对数据进行非线性变换，如log转换和power转换等。数据转换的目的是将一些偏态分布的数据转换成对称分布的数据，并消除异方差性的影响，以满足一些线性分析技术的要求。不同的预处理方法会对统计分析结果产生不同的影响（见表1），在实际应用中，我们应该根据具体的研究目的、数据类型以及要选用的统计分析方法综合考虑，选择适当的预处理方式。例如，Robert A.van den Berg等（2006）通过实际代谢组学数据的分析发现，选用不同预处理方法在很大程度上影响着主成分分析（PCA）的结果，自动尺度化（autoscaling）和全距尺度化（range scaling）在对代谢组学数据进行探索性分析时表现更优，其PCA分析后的结果在生物学上能够得到更合理的解释［11］。

表1 常见的数据预处理方法

单变量分析方法

单变量分析方法简便、直观和容易理解，在代谢组学研究中通常用来快速考察各个代谢物在不同类别之间的差异。代谢组学数据在一般情况下难以满足参数检验的条件，使用较多的是非参数检验的方法，如W ilcoxon秩和检验或Kruskal-Wallis检验，t＇检验也是一种比较好的统计检验方法。

由于代谢组学数据具有高维的特点，所以在进行单变量分析时，会面临多重假设检验的问题。如果我们不对每次假设检验的检验水准α进行校正，则总体犯一类错误的概率会明显增加。一种解决方法是采用Bonferion校正，即用原检验水准除以假设检验的次数m作为每次假设检验新的检验水准（α/m）。由于Bonferion校正的方法过于保守，会明显降低检验效能，所以在实际中更为流行的一种做法是使用阳性发现错误率（false discovery rate，FDR）。这种方法可用于估计多重假设检验的阳性结果中，可能包含多少假阳性结果。FDR方法不仅能够将假阳性的比例控制在规定的范围内，而且较之传统的方法在检验效能上也得到显著的提高［12］。实际中也可以使用局部FDR（用fdr表示），其定义为某一次检验差异显著时，其结果为假阳性的概率。局部FDR的使用，使得我们能够估计出任意变量为假阳性的概率，通常情况下有FDR≤fdr［13］。

除了进行传统的单变量假设检验分析，代谢组学分析中通常也计算代谢物浓度在两组间的改变倍数值（fold change），如计算某个代谢物浓度在两组中的均值之比，判断该代谢物在两组之间的高低表达。计算ROC曲线下面积（AUC）也是一种经常使用的方法［14］。

多变量分析

代谢组学产生的是高维的数据，单变量分析不能揭示变量间复杂的相互作用关系，因此多变量统计分析在代谢组学数据分析中具有重要的作用。总体来说，代谢组学数据多变量统计分析方法大致可以分为两类：一类为非监督的学习方法，即在不给定样本标签的情况下对训练样本进行学习，如PCA、非线性映射（NLM）等；另一类为有监督的学习方法，即在给定样本标签的情况下对训练样本进行学习，如偏最小二乘判别分析（PLS-DA）、基于正交信号校正的偏最小二乘判别分析（OPLS-DA）、人工神经网络（ANN）、支持向量机（SVM）等。其中，PCA、PLS-DA和OPLS-DA是目前代谢组学领域中使用最为普遍的多变量统计分析方法。

PCA是从原始变量之间的相互关系入手，根据变异最大化的原则将其线性变换到几个独立的综合指标上（即主成分），取2～3个主成分作图，直观地描述不同组别之间的代谢模式差别和聚类结果，并通过载荷图寻找对组间分类有贡献的原始变量作为生物标志物。通常情况下，由于代谢组学数据具有高维、小样本的特性，同时有噪声变量的干扰，PCA的分类结果往往不够理想。尽管如此，PCA作为代谢组学数据的预分析和质量控制步骤，通常用于观察是否具有组间分类趋势和数据离群点［15］。在组间分类趋势明显时，说明其中一定有能够分类的标志物。PCA还可以用于分析质控样品是否聚集在一起，如果很分散或具有一定的变化趋势，则说明检测质量存在一定的问题。Zhang Zhiyu等（2010）通过PCA成功区分了骨肉瘤患者和正常人，并发现良性骨肿瘤患者中有两例是异常值［16］。Kishore K.Pasikanti等（2009）利用PCA对尿液膀胱癌代谢组学数据进行分析后观察到质控样品在PCA得分图上紧密聚集，从而验证了仪器检测的稳定性和代谢组学数据的可靠性［17］。

PLS-DA是目前代谢组学数据分析中最常使用的一种分类方法，它在降维的同时结合了回归模型，并利用一定的判别阈值对回归结果进行判别分析。Zhang Tao等（2013）运用PLS-DA技术分析尿液卵巢癌代谢组学数据，成功将卵巢癌患者和良性卵巢肿瘤患者以及子宫肌瘤患者相互鉴别，并鉴定出组氨酸、色氨酸、核苷酸等多种具有判别能力的卵巢癌生物标志物［18］。PLS的思想是，通过最大化自变量数据和应变量数据集之间的协方差来构建正交得分向量（潜变量或主成分），从而拟合自变量数据和应变量数据之间的线性关系［19］。PLS的降维方法与PCA的不同之处在于PLS既分解自变量X矩阵也分解应变量Y矩阵，并在分解时利用其协方差信息，从而使降维效果较PCA能够更高效地提取组间变异信息［20］。当因变量Y为二分类情况下，通常一类编码为1，另一类编码为0或-1；当因变量Y为多分类时，则需将其化为哑变量。通常，评价PLS-DA模型拟合效果使用R2X、R2Y和Q2Y这三个指标，这些指标越接近1表示PLS-DA模型拟合数据效果越好。其中，R2X和R2Y分别表示PLSDA分类模型所能够解释X和Y矩阵信息的百分比，Q2Y则为通过交叉验证计算得出，用以评价PLS-DA模型的预测能力，Q2Y越大代表模型预测效果较好。实际中，PLS-DA得分图常用来直观地展示模型的分类效果，图中两组样品分离程度越大，说明分类效果越显著。代谢组学数据分析中另一种常用的方法是OPLS-DA，它是PLS-DA的扩展，即首先使用正交信号校正技术，将X矩阵信息分解成与Y相关和不相关的两类信息，然后过滤掉与分类无关的信息，相关的信息主要集中在第一个预测成分。Johan Trygg等认为该方法可以在不降低模型预测能力的前提下，有效减少模型的复杂性和增强模型的解释能力［21］。与PLSDA模型相同，可以用R2X、R2Y、Q2Y和OPLS-DA得分图来评价模型的分类效果。Carolyn M.Slupsky等（2010）使用OPLS-DA发现卵巢癌患者、乳腺癌患者、正常人这三者之间的尿液代谢轮廓显著不同，从而推断尿液代谢组学可能为癌症的特异性诊断提供重要依据［22］。

由于代谢组学数据具有高维、小样本的特性，使用有监督学习方法进行分析时很容易产生过拟合的现象。为此，需要使用置换检验考察PLS-DA在无差异情况下的建模效果［23］。该方法在固定X矩阵的前提下，随机置换Y分类标签n次，每次随机置换后建立新的PLS-DA模型，并计算相应的R2Y和Q2Y；然后，与真实标签模型得到的结果进行比较，用图形直观表达是否有过拟合现象。

由于样本量的不足，通常采用上述的交叉验证和置换检验方法作为模型验证方法。而实际中，在样本量允许的情况下，最为有效的模型验证方法即将整个数据集严格按照时间顺序划分为内部训练数据和外部测试数据两部分，利用内部训练数据建立模型，再对外部测试数据进行预测，客观地评价模型的有效性和适用性。

生物标志物的筛选

代谢组学分析的最终目标是希望从中筛选出潜在的生物相关标志物，从而探索其中的生物代谢机制，因此需要借助一定的特征筛选方法进行变量筛选。对于高维代谢组学数据的特征筛选，研究的目的是从中找出对样本分类能力最强或较强的一个或若干个变量。特征筛选方法主要分为三类：过滤法、封装法和嵌入法［24］。过滤法主要是采用单变量筛选方法对变量进行筛选，优点是简单而快捷，能够快速的降维，如t＇检验、W ilcoxon秩和检验、SAM等方法。封装法是一种多变量特征筛选策略，通常是以判别模型分类准确性作为优化函数的前向选择、后向选择和浮动搜索特征变量的算法，它通常是按照“节省原则”进行特征筛选，最终模型可能仅保留其中很少部分的重要变量，如遗传算法等。嵌入法的基本思想是将变量选择与分类模型的建立融合在一起，变量的重要性评价依靠特定分类模型的算法实现，在建立模型的同时，可以给出各变量重要性的得分值，如PLS-DA方法的VIP统计量等。为了更加客观、全面地评价每个变量的重要性，代谢组学研究中一般采取将上述方法结合起来的方式进行变量筛选。比较常见的一种策略是先进行单变量分析，再结合多变量模型中变量重要性评分作为筛选标准，如挑选fdr≤0.05和VIP＞1.5的变量作为潜在生物标志物。用筛选的潜在生物标志物对外部测试数据集进行预测，评价其预测效果。最后，可以通过研究生物标志物的生物学功能和代谢通路，分析不同生物标志物之间的相互作用和关系，从而为探索生物代谢机制提供重要线索和信息。Yang Jinglei等（2013）即在代谢组学分析中使用fdr≤0.2和VIP＞1.5的双重标准来筛选精神分裂症的特异生物标志物，所筛选出的差异代谢物其AUC在训练数据中达94.5%，外部测试数据中达0.895［25］。

展 望

由于代谢组学数据变量多、关系复杂的特性，数据分析任务极为艰巨。目前常用的统计学方法在一定程度上为进行代谢组学数据分析提供了有效的工具，但仍然存在诸多不足。如在代谢组学研究中，生物样品之间的变异性往往较大，目前最流行的PLS-DA或OPLS-DA数据分析方法在差异小、噪声大时，模型使用效果不够理想。另外，PLS-DA和OPLS-DA均是基于线性回归的方法，但是代谢组学数据通常不是简单的线性关系，因此，PLS-DA和OPLS-DA模型拟合数据的结果可能会不够好。基于这些问题，一些学者开始尝试将一些新的高维数据分析方法和思想应用于代谢组学数据分析中，如Lin Xiaohui等（2011）提出一种将支持向量机、随机森林和遗传算法结合起来进行变量筛选的分析思路，通过比较证实其较单个分析方法能够发掘出更多的信息，尤其适合分析复杂生物数据［26］；Elon Correa和Royston Goodacre（2011）提出了一种新型的遗传算法—贝叶斯网络方法（GA-BN），这种方法在有效筛选变量并提高分类效果的同时，还能研究变量间的相互作用和关系［27］。毫无疑问，这些新方法的提出将会为代谢组学数据分析提供新的思路和契机。随着各种代谢组学检测仪器的快速发展，更有效的代谢组学数据分析技术亟待开发，值得更多的生物统计学者关注和研究。

1.Dunn WB，Ellis DI.Metabolom ics：Current analytical platforms and methodologies.TrAC Trends in Analytical Chem istry，2005，24（4）：285-294.

2.许国旺，路鑫，杨胜利.代谢组学研究进展.中国医学科学院学报，2007，29（6）：701-711.

3.Spratlin JL，Serkova NJ，Eckhardt SG.Clinical applications of metabolomics in oncology：a review.Clin Cancer Res，2009，15（2）：431-440.

4.W ishart DS.Applications ofmetabolom ics in drug discovery and development.Drugs R D，2008，9（5）：307-322.

5.Taylor J，King RD，Altmann T，et al.Application of metabolom ics to plant genotype discrim ination using statistics and machine learning. Bioinformatics，2002，18（2）：241-248.

6.Nicholson JK，Lindon JC，Holmes E.＇Metabonom ics′：understanding the metabolic responses of living systems to pathophysiological stimuli via multivariate statistical analysis of biological NMR spectroscopic data. Xenobiotica，1999，29（11）：1181-1189.

7.Smith CA，Want EJ，O′Maille G，etal.XCMS：processingmass spectrometry data formetabolite profiling using nonlinear peak alignment，matching，and identification.Analytical Chemistry，2006，78（3）：779-787.

8.Sima C，Dougherty ER.What should be expected from feature selection in small-sample settings.Bioinformatics，2006，22（19）：2430-2436.

9.Dunn WB，Ellis DI.Metabolom ics：Current analytical platforms and methodologies.Trac-Trend Anal Chem，2005，24（4）：285-294.

10.Goodacre R，Broadhurst D，Sm ilde A，et al.Proposed m inimum reporting standards for data analysis inmetabolom ics.Metabolom ics，2007，3（3）：231-241.

11.Van den Berg RA，Hoefsloot HCJ，Westerhuis JA，etal.Centering，scaling，and transformations：improving the biological information content ofmetabolom ics data.BMC Genom ics，2006，7：142-156.

12.Benjam ini Y，Hochberg Y.Controlling the false discovery rate：a practical and powerful approach to multiple testing.JR Statist Soc B，1995，57（1）：289-300.

13.刘晋，张涛，李康.多重假设检验中FDR的控制与估计方法.中国卫生统计，2012，29（2）：305-308.

14.Broadhurst DI，Kella DB.Statistical strategies for avoiding false discoveries in metabolom ics and related experiments.Metabolom ics，2006，2（4）：171-196.

15.Trygg J，Holmes E，Lundstedt T.Chemometrics inmetabonom ics.JProteome Res，2007，6（2）：469-479.

16.Zhang Z，Qiu Y，Hua Y，etal.Serum and urinarymetabonom ic study of human osteosarcoma.JProteome Res，2010，9（9）：4861-4868.

17.Pasikanti KK，Esuvaranathan K，Ho PC，et al.Noninvasive urinary metabonomic diagnosis of human bladder cancer.J Proteome Res，2009，9（6）：2988-2995.

18.Zhang T，Wu XY，Ke CF，et al.Identification of Potential Biomarkers for Ovarian Cancer by Urinary Metabolom ic Profiling.JProteome Res，2013，12（1）：505-516.

19.蒋红卫，夏结来.偏最小二乘回归及其应用.第四军医大学学报，2003，24（3）：280-283.

20.Boulesteix AL，Strimmer K.Partial least squares：a versatile tool for the analysis of high-dimensional genom ic data.Brief Bioinform.2007，8（1）：32-44.

21.Trygg J，Wold S.Orthogonal projections to latent structures（O-PLS）. Journal of Chemometrics，2002，16（3）：119-128.

22.Slupsky CM，Steed H，Wells TH，et al.Urine metabolite analysis offers potential early diagnosis of ovarian and breast cancers.Clin Cancer Res，2010，16（23）：5835-5841.

23.Westerhuis JA，Hoefsloot HCJ，Sm it S，et al.Assessment of PLSDA cross validation.Metabolom ics，2008，4（1）：81-89.

24.Yvan S，Iaki I，Pedro L.A review of feature selection techniques in bioinformatics.Bioinformatics，2007，23（13）：273-281.

25.Yang J，Chen T，Sun L，et al.Potentialmetabolite markers of schizophrenia.Molecular Psychiatry，2013，18（1）：67-78.

26.Lin XH，Wang QC，Yin PY，etal.A method for handlingmetabonomics data from liquid chromatography/mass spectrometry：combinational use of support vector machine recursive feature elimination，genetic algorithm and random forest for feature selection.Metabolomics，2011，7（4）：549-558.

27.Correa E，Goodacre R.A genetic algorithm-Bayesian network approach for the a nalysis ofmetabolom ics and spectroscopic data：application to the rapid identification of Bacillus spores and classification of Bacillus species.BMC Bioinformatics，2011，12（1）：33-49.

（责任编辑：郭海强）

*国家自然科学基金资助（81172767）；高等学校博士学科专项基金（20122307110004）

1哈尔滨医科大学卫生统计学教研室（150081）

2山东大学卫生统计学教研室

Δ通信作者：李康，E-mail：likang＠ems.hrbmu.edu.cn