代謝組學數據分析的統計學方法*

柯朝甫張 濤武曉巖李 康Δ

代謝組學數據分析的統計學方法*

柯朝甫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