基于色譜-質譜平臺的代謝組學數據預處理方法*

哈爾濱醫科大學衛生統計學教研室(150081)

孫 琳 張秋菊 王文佶 曲思楊 謝 彪 高 兵 劉美娜△

?

·方法介紹·

基于色譜-質譜平臺的代謝組學數據預處理方法*

哈爾濱醫科大學衛生統計學教研室(150081)

孫 琳 張秋菊 王文佶 曲思楊 謝 彪 高 兵 劉美娜△

代謝組學的概念自20世紀90年代被正式提出[1]，已被廣泛應用于醫學研究領域，其一般研究流程包括樣本采集、樣本檢測、數據預處理、數據分析和生物學解釋等。常用的樣本檢測技術有核磁共振(nuclear magnetic resonance，NMR)和高分辨率色譜-質譜聯用技術[2]，本文所述方法針對后者。經色譜-質譜聯用平臺檢測的數據具有以下特點：高維度，小樣本，變量間高度相關，高噪聲，高缺失，以及高度變異性。基于以上數據特征，在代謝組學數據分析之前需要對數據進行預處理[3]，以消除或減小數據中高噪聲、高缺失和高度變異性對統計分析結果的干擾。數據預處理是數據分析的前提，有利于統計分析過程的多變量模型構建；若數據未經過充分預處理而直接用于統計分析，可能會掩蓋數據中某些真實特征(如組間差異)，使研究結果模糊化。數據預處理包括缺失數據處理，數據標準化，以及數據的中心化、標度化和轉換等內容。本文將全面介紹基于色譜-質譜聯用平臺的代謝組學數據預處理方法，并進行各方法比較，為研究者選擇合適的數據預處理方法提供思路。

數據預處理的重要性

代謝組學研究經常涉及處于不同實驗條件下兩組樣本的分類判別(例如處理組和非處理組、病例組和對照組)，目的是建立某種疾病或危險因素的預測模型，找到相關的差異表達生物標志物，輔助疾病診斷、預防或治療[4]。代謝組學數據中存在多種變異，其中，由于某種特定誘導因素(實驗干預、疾病等)導致的代謝物濃度差異被稱為誘導變異，這是研究者感興趣的變異，也是數據分析所針對的變異；此外，實驗過程中其他因素引起的變異也會影響檢測的代謝物峰強度，從而影響研究對象分類，主要包括非誘導變異和技術平臺變異。

非誘導變異[5]包括：(1)量級差異：指機體內某個單分子物質的平均濃度遠小于高濃度化合物(例如ATP)的平均濃度，從生物學角度來說，相對于低濃度的代謝物，高濃度代謝物在機體生理病理過程中未必發揮更重要的作用；(2)倍數差異：指誘導因素對不同通路的代謝物的影響程度不一致，一般來說，中心代謝途徑物質的濃度較為穩定，而次級代謝途徑產生的物質的濃度易受環境因素影響；(3)生物固有差異：指在相同的實驗條件下，某些代謝物的濃度在同組別的個體間也展現出較大幅度波動，這種個體差異可能會導致錯誤的判別結果。技術平臺變異來源于樣本預處理過程和色譜-質譜檢測平臺，例如預處理時樣本間細微的體積差異、質譜離子源的波動和色譜柱效能的變化等。檢測平臺變異導致測量誤差，具體表現為數據中的噪聲，在統計分析時，通常假定代謝物噪聲強度呈正態分布，在實際研究中這種假設往往不成立，偏態分布的噪聲造成數據的異方差性，影響統計分析結果。

數據預處理的目的就是減少諸多不利因素對多變量模型建立和生物標志物篩選的干擾，以保證研究結果的準確性。代謝特征數據由色譜-質譜平臺輸出后，要經過原始數據的預加工和清潔數據的預處理，才能進入數據分析階段。

清潔數據

在數據預處理前，首先需要將原始數據轉換為清潔數據。生物樣本經色譜-質譜聯用平臺檢測得到由質荷比、保留時間和代謝物相對峰強度信息構成的三維原始數據，這種數據雜亂且不匹配，不能直接反映研究對象的代謝特征。首先要對原始數據進行預加工，包括保留時間對齊、峰提取和峰匹配等步驟，把原始數據轉換成匹配的色譜-質譜-代謝物數據。預加工后的數據矩陣由幾十到上百個樣本(觀測)所對應的上千個代謝物(變量)的質譜峰強度構成，稱為清潔數據。實際研究中，預處理和統計分析階段的原始數據指的是清潔數據而非儀器導出的原始數據。圖1展示了整個數據預處理的流程。

常用的代謝組學數據預處理方法

1．缺失值處理 色譜-質譜聯用平臺輸出數據中存在大量“0”值，我們把這些“0”值標記為缺失數據。產生缺失值的原因有三方面[6]：代謝物只存在于某些樣本，在其他樣本中濃度為0；代謝物存在于樣本中，但其濃度低于檢測儀器設定的檢出限；代謝物在樣本中的濃度高于檢出限而被色譜檢出，但沒有得到正確的質譜峰匹配，在質譜中無法檢出。缺失值的產生會對深入數據挖掘結果產生較大偏倚，降低統計分析的效率，因此需要對缺失數據進行刪減和填補。

圖1 代謝組學數據預處理流程圖

(1)缺失值刪減 一般遵守“80規則”[6]，即：如果某一代謝物在一個類別超過80%的樣本中的檢測強度不為0，則該變量予以保留；反之剔除。通過“80規則”，由于檢測儀器導致的數據缺失被移除。

(2)缺失值填補 在初步缺失值刪減后，數據中仍存在一部分“0”值，需要對這部分數據進行填補。目前常用的缺失值填補方法包括列替代法[7-8]、K鄰近法[8-9]、貝葉斯主成分分析法[7]和多重填補法[10-12]等。

列替代法 由于研究對象的同質性，代謝物在同類別研究對象的不同個體間的代謝行為高度相似。因此，某種代謝物在特定實驗條件下的缺失值，可用缺失數據所在列的其他樣本的數據估計缺失的代謝物強度，包括最小值替代法、均值替代法和中位數替代法。列替代法簡單易操作，但是并沒有考慮到數據中各變量的相關性，因此估計的準確度較低。

K鄰近法 首先計算含有缺失值的代謝物和所有其他代謝物間的歐氏距離，選取與缺失變量歐氏距離最近的K個變量的加權平均峰強度值填補缺失數據， K選擇10至20較為合理。不同于列替代法，K鄰近填補法對每個缺失值給出不同估計值，更接近真實數據，而且考慮了變量間的相關性。

貝葉斯主成分分析法(Bayesian principal components analysis，BPCA) 此種填補方法是基于主成分回歸、貝葉斯估計和EM算法的三層算法。先用主成分分析的結果建立完整數據的主成分回歸模型，用回歸模型的預測值作為缺失數據的初步估計值，建立貝葉斯估計的先驗分布，假定殘差和代謝物在主成分上的投影都是正態獨立的變量，分布參數未知，由EM算法迭代估計上述參數直至收斂，用最終的主成分回歸模型預測值作為缺失數據的估計值。BPCA填補的優勢在于利用了數據的全部結構，對樣本量較大的隨機缺失處理效果較好[28-29]。

多重填補法 與單一填補不同，多重填補用某個缺失值的可能取值的集合進行填補，重復多次，產生多個完整數據集；隨后對多個數據集進行統計分析，得到每個缺失數據的平均水平和變異水平；最后，整合來自各填補數據集的結果，從而得到統計推斷的結果。該方法考慮到了缺失數據的不確定性對統計推斷結果的影響。

以上介紹的方法側重點不同，實際研究中，判斷是否需要進行缺失值填補以及選擇何種填補方法需要研究者依據數據特征決定。

2．數據標準化(表1) 標準化針對代謝組學實驗中的系統誤差[13]，包括樣本采集過程中的差異(如尿液樣本采集時間不同)、個體差異(如由于飲水量不同而導致的尿液濃度差異)或大規模研究中的批次差異。與其他數據預處理方法相比，標準化方法更為復雜 ，這里簡述幾種常用的標準化方法。

(1)總強度標準化 假定數據中所有樣本的全部特征的峰強度總和相等，即所有代謝物的總強度在樣本之間不發生變化。每個觀測中所有元素的標準化因子相同，將原始數據除以所有特征的峰強度之和得到新數據[14]。這是代謝組學中最常用的標準化方法，Waters公司的軟件MakerLynx中默認標準化方法就是總強度標準化。

(2)Cyclic Loess標準化 其基本思想是通過曲線擬合的方法，調整兩個樣本間的所有特征的對數峰強度值差值，使之近似為0[15-16]。Cyclic Loess對不同類別樣本數據分別進行標準化，可以有效地減少組間個體變異和技術平臺變異。

(3)Quantile標準化 是基因組學中常用的標準化方法[18]，目的是使全部樣本的數據具有相同分布[17]，當利用Q-Q圖對數據進行可視化時，所有點近似分布于對角線上。通過使各樣本數據分位數相等，達到所有觀測代謝物數據的一致分布。

3.數據的中心化、標度化和轉換(表2-3)

對于非誘導變異和技術變異，除上述標準化方法外，還可以采取中心化和標度化處理。轉換則是對由于噪聲導致的數據異方差性進行處理。通過這一階段的預處理，數據中來源于誘導因素的代謝物濃度變異得以突顯。

表1 標準化方法總結

(1)數據中心化

具體方法是將原數據減去每個列變量的均值而得到新數據，中心化后，新數據圍繞0上下波動，而不再圍繞代謝物峰強度均值波動。通過中心化，高濃度代謝物和低濃度代謝物間的濃度差異得到適當的調整，突出數據中波動部分的重要性[19-20]。中心化是標度化和轉換的基礎，以下提到的兩類方法都要與中心化相結合。

(2)數據標度化

Auto scaling 自標度化，又稱為單位標度化或單位方差標度化。以標準差作為標度因子，用中心化后的數據除以列變量的標準差而得到新數據[20]。自標度化后，所有代謝物強度的標準差均為1，相當于把所有的變量置于同等重要的水平，消除了由于代謝物的絕對含量差異引起的高濃度物質對低濃度物質的數據掩蓋。需要注意的是，自動標度化給所有代謝物以相同權重，使低含量代謝物更可比，但是，這些物質的測量誤差同時被放大 ，可能會導致錯誤的數據分析結果。

Range scaling 極差標度化，采用生物學范圍作為標度因子[21]，這里的生物學范圍指的是某種代謝物強度的最大值與最小值之差，即列變量的極差。標準差涵蓋了所有研究對象的信息，而極差的計算只用到兩個樣本的信息， 因此極差標度化對于數據中的離群值很敏感。為了提高方法的魯棒性，可以采用其他更穩健的生物學范圍計算方法。

Pareto scaling 該方法和自動標度化類似，只是標度因子由標準差變為標準差的平方根[22]。Pareto標度化在弱化高含量代謝物重要性的同時又保持了數據的完整性，與自動標度化方法相比，新數據更接近原始測量值。與在組別間變化不顯著的代謝物相比，變異較大的物質對此方法更加敏感。

Vast scaling 大規模標度化，又稱為變量穩定性標度化，是自標度化的擴展[23]。大規模標度化重點關注那些在樣本間變異較小，濃度較穩定的代謝物，采用標準差和變異系數作為標度因子。變異系數是列變量標準差和均數的比值，引入變異系數作為標度因子，提高了低相對標準偏差變量的重要性，降低了高相對標準偏差變量的重要性。大規模標度化既可用于無監督模式識別，也可用于有監督模式識別，當進行有監督模式識別時，用各組的信息分別求組內變異系數作為標度因子。

Level scaling 水平標度化，使用的標度因子是代謝物的平均強度[5]，標度化后的新數據是某種代謝物強度相對于其平均強度變化的百分比，即相對變化值。水平標度化適用于研究特定的某些相對變化大的生物學反應(如應激反應)和高濃度生物標志物的發現。

(3)數據轉換

轉換是非線性的數據預處理方法，主要用于校正數據中的噪聲結構[24]，通過數據轉換，將乘性噪聲轉換為加性噪聲，使偏態分布的數據轉換為更加對稱的分布。

經UPLC-MS檢測輸出的峰強度包括代謝物的信號強度以及噪聲強度這兩部分，其中噪聲具有乘性結構或加性結構[25-26]，可以用以下模型表示：

xij=βij+nij·sijeηij

xij是檢測到的第i種代謝物在第j個樣本中的峰強度，sij是預期的峰強度，βij表示隨機背景噪聲(電子噪聲)，ηij表示乘性隨機噪聲(例如樣本預處理過程出現的變異，離子源的波動或樣品導入設備的波動)，ηij表示樣本j的標準化因子。加性噪聲是基線的隨機波動，與代謝物的信號強度無關，來自于檢測儀器的電子噪聲，由于交叉進樣，加性噪聲均勻出現在各組別的全部樣本的檢測過程中，因此可忽略不計。與此相反，乘性噪聲隨代謝物信號強度增加而增強，通常噪聲強度與信號強度成比例。由于乘性噪聲的存在，進行不同樣本的多次測量時，高濃度的代謝物會展現出更大的變異，造成某些低濃度代謝物的信號可能會被高濃度代謝物的強噪聲所掩蓋。

表2 數據中心化和數據標度化方法的總結

對數轉換 是代謝組學常用的數據轉換方法。如果數據中各變量的相對標準偏差是一個常量，則使用對數轉換能完全消除乘性噪聲對代謝物峰強度的累加作用[24]，在實際研究中，這種情況非常罕見。對數轉換無法處理數據中的零值，因此對數轉換前先要對數據進行缺失值填補。對數轉換的另一個缺陷是不能很好地處理那些峰強度相對標準偏差很大的代謝物，這些物質通常有較低的相對濃度，變異相對于均數更加突出，當峰強度趨于0時，對數轉換后的數值接近負無窮。

平方根轉換 對于缺失值較多或低濃度代謝物較多的數據，通常選擇平方根轉換[27]來校正噪聲強度。平方根轉換后，高濃度代謝物的方差明顯減小，使得濃度不同的代謝物的方差近似相等。平方根轉換只是縮小了乘性噪聲強度，并沒有將其轉變為加性噪聲。

表3 數據轉換方法總結

小 結

在代謝組學研究過程中，數據預處理會對統計分析結果產生很大影響，因此數據分析前的預處理過程是必不可少的。預處理方式多樣，關于預處理方法的比較研究較少 。本文介紹了多種數據缺失值填補法，標準化法以及標度化、轉換法。

從理論上看，BPCA填補法和多重填補法對缺失值的估計利用了數據集中全部的信息，填補效果更好。Hrydziuszko等[30]比較各種缺失值填補方法(包括基于單變量方法和基于多變量的方法)處理代謝組學數據的效果，發現K鄰近法的填補效果最優。

Bedilu[31]等比較不同標準化方法對于大規模代謝組學數據的影響，發現曲線擬合的Cyclic Loess標準化法能最有效地移除系統誤差。Kévin[32]等通過對實驗數據的分析發現，在幾種標準化方法中，Cyclic Loess效果較好，線性的標準化法效果較差。Pirtr S[33]等比較不同標度化方法基于主成分-判別分析、隨機森林和K最鄰近分類三種多變量分析方法對于統計分析結果的影響，發現自動標度化和極差標度化優于其他方法，Robert等[5]也得到了類似的結論。

數據轉換對大數值的影響相較于小數值更大，縮小了高濃度代謝物和低濃度代謝物濃度之間的量級差異，起到了類似標度化的偽標度化作用。大規模大樣本量的代謝組學研究涉及到多批次數據的合并，數據中的噪聲結構更為復雜，可以考慮標度化和轉換兩種方法的結合運用。在高維組學數據中，數據轉換和標度化之間的相互影響很復雜，應用何種方法還要依據實際情況而定。

預處理方法選擇不僅取決于代謝數據的特點，還與之后選用的數據分析方法有關。例如，聚類方法關注于特征間相似點(或不同點)的分析，而主成分分析則是試圖用盡可能少的成分解釋數據中大部分變異(降維)；數據預處理可能會改善聚類方法的結果，卻使主成分分析的結果模糊化 。值得注意的是，上述預處理方法比較研究的結論都基于研究者在特定條件下模擬得到的數據集和特定的實驗數據集的分析而得出的，在外推到其他數據集時要慎重。選擇數據預處理方法時要綜合考慮研究目的、數據結構特征和擬采用的統計分析方法，再決定預處理的策略和具體方法。

[1]Oliver SG,Winson MK,Kell DB,et al.Systematic functional analysis of the yeast genome.Trends in biotechnology,1998,16(9):373-378.

[2]Salek RM,Steinbeck C,Viant MR,et al.The role of reporting standards for metabolite annotation and identification in metabolomic studies.GigaScience,2013,2(1):1.

[3]Goodacre R,Broadhurst D,Smilde AK,et al.Proposed minimum reporting standards for data analysis in metabolomics.Metabolomics,2007,3(3):231-241.

[4]Fiehn O.Metabolomics-the link between genotypes and phenotypes.Plant molecular biology,2002,48(1-2):155-171.

[5]van den Berg RA,Hoefsloot HCJ,Westerhuis JA,et al.Centering,scaling,and transformations:improving the biological information content of metabolomics data.BMC genomics,2006,7(1):1.

[6]Bijlsma S,Bobeldijk I,Verheij ER,et al.Large-scale human metabolomics studies:a strategy for data(pre-)processing and validation.Analytical chemistry,2006,78(2):567-574.

[7]Xia J,Psychogios N,Young N,et al.MetaboAnalyst:a web server for metabolomic data analysis and interpretation.Nucleic acids research,2009,37(suppl 2):W652-W660.

[8]Steuer R,Morgenthal K,Weckwerth W,et al.A gentle guide to the analysis of metabolomic data.Metabolomics:Methods and protocols,2007:105-126.

[9]Troyanskaya O,Cantor M,Sherlock G,et al.Missing value estimation methods for DNA microarrays.Bioinformatics,2001,17(6):520-525.

[10]Rubin DB.Multiple imputations in sample surveys-a phenomenological Bayesian approach to nonresponse//Proceedings of the survey research methods section of the American Statistical Association.American Statistical Association,1978,1:20-34.

[11]Schafer JL.Multiple imputation:a primer.Statistical methods in medical research,1999,8(1):3-15.

[12]Buuren S,Groothuis-Oudshoorn K.mice:Multivariate imputation by chained equations in R.Journal of statistical software,2011,45(3)：.

[13]Kultima K,Nilsson A,Scholz B,et al.Development and evaluation of normalization methods for label-free relative quantification of endogenous peptides.Molecular & Cellular Proteomics,2009,8(10):2285-2295.

[14]Veselkov KA,Vingara LK,Masson P,et al.Optimized preprocessing of ultra-performance liquid chromatography/mass spectrometry urinary metabolic profiles for improved information recovery.Analytical chemistry,2011,83(15):5864-5872.

[15]Cleveland WS,Devlin SJ.Locally weighted regression:an approach to regression analysis by local fitting.Journal of the American Statistical Association,1988,83(403):596-610.

[16]Dudoit S,Yang YH,Callow MJ,et al.Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments.Statistica sinica,2002:111-139.

[17]Bolstad BM,Irizarry RA,?strand M,et al.A comparison of normalization methods for high density oligonucleotide array data based on variance and bias.Bioinformatics,2003,19(2):185-193.

[18]Hansen KD,Irizarry RA,Zhijin WU.Removing technical variability in RNA-seq data using conditional quantile normalization.Biostatistics,2012,13(2):204-216.

[19]Bro R,Smilde AK.Centering and scaling in component analysis.Journal of Chemometrics,2003,17(1):16-33.

[20]Jackson JE.A user's guide to principal components.1991.

[21]Smilde AK,van der Werf MJ,Bijlsma S,et al.Fusion of mass spectrometry-based metabolomics data.Analytical chemistry,2005,77(20):6729-6736.

[22]Eriksson L.Introduction to multi-and megavariate data analysis using projection methods(PCA & PLS).Umetrics AB,1999.

[23]Keun HC,Ebbels TMD,Antti H,et al.Improved analysis of multivariate data by variable stability scaling:application to NMR-based metabolic profiling.Analytica chimica acta,2003,490(1):265-276.

[24]Kvalheim OM,Brakstad F,Liang Y.Preprocessing of analytical profiles in the presence of homoscedastic or heteroscedastic noise.Analytical Chemistry,1994,66(1):43-51.

[25]Anderle M,Roy S,Lin H,et al.Quantifying reproducibility for differential proteomics:noise analysis for protein liquid chromatography-mass spectrometry of human serum.Bioinformatics,2004,20(18):3575-3582.

[26]Durbin BP,Hardin JS,Hawkins DM,et al.A variance-stabilizing transformation for gene-expression microarray data.Bioinformatics,2002,18(suppl 1):S105-S110.

[27]Sokal RR,Rohlf FJ.Assumptions of analysis of variance.Biometry:The principles and practice of statistics in biological research,1995:392-450.

[28]Albrecht D,Kniemeyer O,Brakhage AA,et al.Missing values in gel-based proteomics.Proteomics,2010,10(6):1202-1211.

[29]Oba S,Sato M,Takemasa I,et al.A Bayesian missing value estimation method for gene expression profile data.Bioinformatics,2003,19(16):2088-2096.

[30]Hrydziuszko O,Viant MR.Missing values in mass spectrometry based metabolomics:an undervalued step in the data processing pipeline.Metabolomics,2012,8(1):161-174.

[31]Ejigu BA,Valkenborg D,Baggerman G,et al.Evaluation of normalization methods to pave the way towards large-scale LC-MS-based metabolomics profiling experiments.Omics:a journal of integrative biology,2013,17(9):473-485.

[32]Contrepois K,Jiang L,Snyder M.Optimized analytical procedures for the untargeted metabolomic profiling of human urine and plasma by combining hydrophilic interaction(HILIC)and reverse-phase liquid chromatography(RPLC)-Mass spectrometry.Molecular & Cellular Proteomics,2015,14(6):1684-1695.

[33]Gromski PS,Xu Y,Hollywood KA,et al.The influence of scaling metabolomics data on model classification accuracy.Metabolomics,2015,11(3):684-695.

(責任編輯：鄧 妍)

國家科技支撐計劃(2011BAI09B02)；黑龍江省自然基金重點項目(ZD201314)

△通信作者：劉美娜，E-mail： liumeina369@163.com