動態代謝組學數據分析方法介紹*

哈爾濱醫科大學公共衛生學院衛生統計學教研室（150081） 王文佶 張秋菊 曲思楊 謝彪 高兵 孫琳 劉美娜

動態代謝組學數據分析方法介紹*

哈爾濱醫科大學公共衛生學院衛生統計學教研室（150081） 王文佶 張秋菊 曲思楊 謝彪 高兵 孫琳 劉美娜△

近年來，代謝組學發展迅速并廣泛應用于營養學、毒理學、疾病診斷等各個領域［1］。隨著研究的深入，代謝組學所獲得的數據集因研究設計的不同而日益復雜，產生了有時間間隔的動態代謝組學數據。然而目前分析此類數據的方法十分有限，并且在大多數情況下，這種動態數據所帶來的因時間因素產生的變異直接被忽略。據此，本文回顧了現有的動態代謝組學數據的分析方法并對其研究進展進行介紹。

常用代謝組學數據分析方法

代謝組學數據最大特點是變量數遠遠大于觀測數，且變量之間存在著高度相關性。目前，最常用于代謝組學數據分析的方法有：無監督學習方法的主成分分析（PCA）［2］、有監督學習方法的偏最小二乘判別分析（PLS-DA）［3］、正交偏最小二乘判別分析（OPLSDA）等［4］。這些方法可以提取原始信息的最大變異或在此基礎上的最佳解釋變異，將高維數據映射到低維空間，并給出降維后數據的可視化展示。隨著研究深入，代謝組學不再拘泥于橫斷面研究，包含時間因素的動態代謝組數據被越來越多的研究所獲得，這使得傳統分析方法不再適用［5］。以PCA為例，來說明傳統分析方法不再適用動態代謝組學數據分析的原因。

在動態代謝組學研究中，其數據特點不止是小樣本大變量所帶來的高維性，數據中還存在與時間有關的變異。例如：研究某種干預（藥物，毒物等）隨著時間改變對不同處理組產生的效應是否有差別，并感興趣于找出隨時間改變的差異代謝物。若以矩陣每一行代表代謝組學數據每一個樣本的觀測時間點，列代表代謝產物進行PCA，在不斷進行行間打亂后，會發現原始得分矩陣Z始終等同于打亂后的得分矩陣Zr，原始載荷矩陣P0始終等同于打亂后的載荷矩陣這說明PCA盲于觀測之間因時間因素產生的自相關變異，若生硬地把所有時間點的數據同時進行主成分分析，只是一味地提取原始信息的最大變異進行降維而忽視變異的來源，將導致與時間相關的變異成為混雜因素而混淆處理因素帶來的差異。這就是傳統PCA、PLS-DA、OPLS-DA不能解決動態代謝組學的關鍵。此時，迫切需要一類代謝組學數據的分析方法，可以在降維的同時捕捉到時間趨勢產生的變異，更真實地揭示數據之間內部結構——動態代謝組學數據分析方法。

動態代謝組學數據分析方法

本文主要從降維的角度對動態代謝組學數據分析方法進行歸納介紹，包括以下方法：動態主成分分析、動態偏最小二乘法、方差同步主成分分析以及動態概率主成分分析方法等。從模型的發展來源、基本原理、適用情況及優缺點等方面介紹。

1．動態主成分分析

動態PCA模型包括兩種形式，第一種是由Ku等人1995年提出的動態PCA［7］，第二種由Timmerman等人在2001年提出的滯后PCA［8］，區別在于前者是對X矩陣進行二次變換，后者是對得分矩陣進行二次變換。其基本思想是：每個觀測在每一個時間點所測得的代謝物濃度（位移）不僅由本時間點的測量所決定，還受之前時間點影響。所以動態PCA可以看成是自回歸滑動平均外生（ARMAX）的時間序列模型與PCA模型的結合［9］。定義二次變換矩陣Bl＝［0T×（L－l）｜IT｜0T×l］，l＝0，…，L代表時間滯后，X為（（T+L）×J）的矩陣。下面以L＝2來闡述動態PCA的工作原理。

新矩陣實質上是一個三維矩陣，包含三種變異：不同變量（代謝產物）之間的變異；相同變量不同時間點間的變異；不同時間點以及不同變量間的變異［10］。因此對矩陣進行PCA分析后，得到的是這三種變異的混合得分。

動態PCA的一個局限是：由于把X矩陣分割成了各個部分，降低了在時間方向上樣本的數量，導致隨著時間點增多丟失的信息隨之增加。另一種可選擇的三維矩陣分析方法是PARAFAC或Tucker3模型［10］，這兩類模型沒有對變量進行滯后變換，一定程度上解決了動態PCA因多次滯后變換導致信息丟失的問題。Rubingh等人利用PARAFAC和N-PLS對微生物發酵過程中細胞代謝產物進行動態分析，得到了發酵的兩個不同的階段，識別出隨著時間變化，誘導大腸桿菌產生苯丙氨酸的代謝物［11］。

此外，PCA的擴展模型還包括2004年Jansen等人提出的加權PCA［12］，其基本思想是：根據代謝物同一時間點測量值的變異大小匹配不同的權重，除去了測量誤差的影響，凸顯代謝物在時間上的波動性。2009年Jansen等人在一篇關于植物代謝組學研究中［13］，提出局部PCA模型：先對所有時間點進行整體的PCA分析，再把每個時間點的變量應用局部PCA模型，最后利用交叉適應分析（cross fit analysis）將幾個模型不同時間點聯系起來作為動態分析結果解釋。

2．動態偏最小二乘法

這里介紹三種動態偏最小二乘法：

第一種方法是有限脈沖序列模型，1987年Ljung等人在系統識別應用中提出［14］。該方法首先對X矩陣進行時間滯后變換，然后將滯后的x變量作為一個新矩陣與表型y進行偏最小二乘分析（PLS）。以此模型為基礎，1996年Qin等人在滯后變量x的基礎上同時滯后變量y，再進行建模分析，對該模型進行了擴展［15］。然而，基于代謝組學數據小樣本多變量的特點，這兩種方法使得原本已經龐大的X矩陣，在不斷被滯后變量擴展后變得更為龐大，導致低樣本量帶來的問題進一步加劇，且增大的載荷矩陣讓模型很難解釋，無法挖掘潛在的數據結構［16］。此時PLS雖有著較好的降維能力，但與其帶來的缺點相比，此方法并不可行［16］。

第二種方法于1993年由Kaspar和Ray等人在過程控制中提出［17］。該方法的核心是利用先驗的動態信息定義了一個動態“過濾器”，將濾過的X矩陣與y建立PLS模型。這種方法不存在上述X維度突長的問題，但在代謝組學中，對這樣敏感的“過濾器”的調節是一件艱巨而又復雜的任務。

第三種方法為批次建模（batchmodeling），于1998年由Wold等人提出。此方法最早應用在化工領域來監視批次生產過程［18］。方法的基本思想是把時間當做一個額外的變量參與建模，所研究數據往往是三個維度，各個維度分別代表：觀測對象、變量、樣本的時間點。Antti等人用該方法在小鼠尿液核磁共振組學的研究中［19］，得到肝毒素動物模型中與時間有關的代謝變化。該方法的缺點是所有的研究對象必須有類似的代謝變化和反應速率。此外，2008年Rantalainen等人提出了另一種處理代謝組學時間序列數據的新方法——分段多元建模［20］，該方法保留了正交偏最小二乘的特點，并且可以描述相鄰時間點的差異，但若兩個相鄰時間點代謝物的差異較小，或取樣的頻率較高時，該方法變得不再適用。

3．方差同步主成分分析（ANOVA-simultaneous component analysis，ASCA）

處理動態數據時也可選用方差分析（ANOVA）模型，模型中時間因素可以以定性或定量的方式來估計［21］。定性分析是指時間因素被當成一個定性變量擬合至模型中，從而得到時間效應。定量分析是指時間因素被當成一個定量變量以線性、二次方或曲線的形式擬合到模型中。在有關基因表達的文獻中，我們可以看到時間因素被定性和定量擬合的實例［22－23］。

由于ANOVA只適合單一的生物反應，在面對代謝組學的高通量數據時，忽略了代謝產物之間的高相關性，因此不再適用。關于ANOVA方法的擴展層出不窮，1979年Mardia等人提出多元方差分析模型（multivariate-ANOVA，MANOVA）［24］，但因代謝組學數據協方差矩陣的奇異性使得該方法步履維艱。2004年Smilde等人提出方差同步主成分分析（ANOVA-simultaneous component analysis，ASCA）［25］，把數據集分解成效應矩陣同時進行降維分析。以基于實驗設計、考慮時間因素的代謝數據為例，ASCA首先將總變異分解成三個來源：時間因素，時間因素與處理因素的交互效應，以及不同個體間的變異。分別對這三部分矩陣進行主成分分析，得到隨著時間的改變各部分效應的得分變化趨勢圖，從而給出合理的生物學解釋。Smilde等人在一篇關于代謝組學干預研究中利用ASCA模型檢驗了維生素C對幾內亞豬骨關節炎進展中的作用［25］。此外，ASCA模型已成功應用到心理測量學［26］、蛋白質組學［27］、以及化學計量學［28］。

ASCA模型還可以與平行因子分析（PARAFAC）相結合，即PARAFASCA模型，與ASCA相比，PARAFASCA模型在描述因實驗設計帶來的交互效應時更加簡潔易懂［29］。此外，ASCA沒有限定時間因素必須為線性，Sm ilde等人用其擴展的方法捕捉非線性時間的變化［30］。需要注意的是我們這里說ASCA是一種動態的方法，其前提是時間因素被定量地擬合到模型中［6］。

4．動態概率主成分分析（dynamic probabilistic principal components analysis，DPPCA）

在代謝組學的數據分析中，不論是適合橫斷面分析的傳統PCA，還是基于代謝組學實驗設計，考慮時間因素所擴展的動態PCA模型、加權PCA、局部PCA、ASCA模型，從方法的角度來說，它們都不是一個相應的生成概率模型，因此很難評價在擬合模型估計時的不確定性；從應用的角度來說，縱向代謝組學研究更關注的是處理因素對代謝物在時間效應上的影響，找出真正隨時間變化有顯著意義的代謝產物，從而挖掘出更多的生物學信息，給出更科學合理的生物學解釋，而上述方法都在模型的解釋性上略顯貧乏。據此，2013年Nyamundanda等人提出動態概率主成分分析（DPPCA），并結合線性混合模型（linearm ixed model，LMM）［31］完美解決了上述問題［1］。

DPPCA模型實質上是概率主成分分析（probabilistic principal components analysis，PPCA）［32］與隨機波動模型（stochastic volatility，SV）［33］的結合。前者主要用于數據降維和分析數據中變量的變異，后者假設變量在各時間點間存在一階自回歸過程，分析由于重復測量而產生的時間變異。PPCA由Tipping和Bishop在1999年基于高斯潛變量模型提出［32］，模型表達式為：xi＝W ui+εi（W是載荷矩陣，ui是潛變量，εi是誤差項）。它將主成分分析理論放在概率框架中進行討論，給出了數據信息在主子空間表達時所對應的概率密度估計。

DPPCA模型的構建原理為：假設PPCA模型中潛變量ui和誤差εi滿足隨機波動模型，觀測在時間點m的p個觀測變量xim的DPPCA模型表達式如下所示：

Wm為載荷矩陣為潛在得分，且Wm和隨著時間而變化。然后分別對DPPCA的潛變量uim和誤差εim建立SV模型，得到個體i（i＝1，…，n），主成分j（j＝1，…，q），時間點m（m＝1，…，M），潛變量uijm的SV模型以及觀測i在時間點m的p維向量的誤差項的SV模型。

縱向代謝組學數據利用DPPCA模型，實現了以下至關重要的幾點：首先在數據可視化上，DPPCA在降維的空間里給出不同處理組的主成分隨時間改變的軌跡。其次，DPPCA可以評價不同處理組的時間效應，通過后驗分布的載荷大小得到不同處理組最具影響力的代謝物清單。最后利用線性混合模型在這些最具影響力的代謝產物中，識別濃度隨時間變化而變化顯著的物質。2013年Nyamundanda等人利用DPPCA模型找到了在給藥組和對照組隨著時間的改變，代謝趨勢相反的代謝產物［32］。

5．其他方法：

除上述降維的方法外，2006年Tai和Speed提出多元貝葉斯時間序列分析法應用于基因芯片的時間序列分析中［34］，Metaboanalyst網站（http：//www．metaboanalyst．ca）基于該方法的改進給出了時間序列代謝組學的實現途徑［35］，模型按照隨著時間改變，組間差異最大的順序得到代謝產物的排名清單。2011年Berker等人構建了一個適合短時間序列的代謝組學數據分析框架［36］，通過平滑樣條混合效應（smoothing splinem ixed effects，SME）模型，使觀測的每個時間序列被視為潛在的隨機分布或平滑曲線被估計，最后用一個函數統計量來檢驗兩組曲線的差別，從而找出隨時間改變不同組間的差異代謝物。此外，聚類算法也被應用于時間序列數據分析中［37－38］，實現通過動態變化來進行特征分組，讀者有興趣可以自行查閱。

小 結

動態代謝組學數據不同于橫斷面研究數據，其數據集的變異來源更加復雜。傳統PCA只能提取數據集的最大變異而無視變異的來源，與時間有關的變異會成為一個混雜因素混淆處理因素的作用，因此選擇一種動態代謝組學數據的分析方法尤為必要。本文主要從降維的角度對動態代謝組學數據分析方法進行回顧與介紹，揭示時間效應對生物代謝的影響。

在選擇一種動態分析模型時應考慮實驗設計所涉及的時間點的個數、時間的間隔、所測得的代謝物的數量。可以根據研究設計和目的不同選用加權PCA，局部PCA批次模型從不同角度解釋時間帶來的變異性。若想檢驗數據不同來源的變異隨時間的變化趨勢可選用ASCA模型；若研究關注隨時間改變不同處理組的差異代謝產物可選用DPPCA和混合線性模型；另外，多時點的動態代謝組學數據慎重選用動態PCA避免因時間點增多丟失的信息逐漸增加。

由于此類方法在降維的同時需要考慮時間因素的變異，使得模型難構建且可能涉及多個模型的結合問題，導致分析方法和分析技術的復雜性增加，限制了該類方法的廣泛應用。一種簡潔而有效的模型亟待開發，值得更多的生物統計學者的研究與關注。一個合適的動態代謝組學數據分析方法可以挖掘出數據潛在的更真實、有價值的信息，對健康、疾病有關的生物節律的探索有著重要意義。

［1］Nyamundanda G，Gorm ley IC，Brennan L．A dynam ic probabilistic principal componentsmodel for the analysis of longitudinalmetabolom ics data．Journal of the Royal Statistical Society：Series C（Applied Statistics），2014，63（5）：763-782．

［2］Bro R，Sm ilde AK．Principal component analysis．Analytical Methods，2014，6（9）：2812-2831．

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

［4］Boccard J，Rutledge DN．A consensus orthogonal partial least squares discrim inant analysis（OPLS-DA）strategy formultiblock Om ics data fusion．Anal Chim Acta，2013，769：30-39．

［5］Sm ilde A．Analysis of High-dimensional Data from Designed Metabolom ics Studies．Metabolic Profiling：Disease and Xenobiotics，2014，21：117．

［6］Sm ilde AK，Westerhuis JA，Hoefsloot HC，et al．Dynam ic metabolom ic data analysis：a tutorial review．Metabolom ics，2010，6（1）：3-17．

［7］Ku W，Storer RH，Georgakis C．Disturbance detection and isolation by dynam ic principal component analysis．Chemometrics and intelligent laboratory systems，1995，30（1）：179-196．

［8］Timmerman ME．Component analysis of multisubject multivariate longitudinal data．Ph D thesis，University of Groningen，2001．

［9］Chen J，Liu K．On-line batch processmonitoring using dynam ic PCA and dynam ic PLS models．Chem ical Engineering Science，2002，57（1）：63-75．

［10］Sm ilde AK，Hendriks MM，Westerhuis JA，et al．Data Processing in Metabolom ics．Metabolom ics in Practice：Successful Strategies to Generate and Analyze Metabolic Data，2013：261-284．

［11］Rubingh CM，Bijlsma S，Jellema RH，et al．Analyzing longitudinal microbialmetabolomics data．J Proteome Res，2009，8（9）：4319-4327．

［12］Jansen JJ，Hoefsloot HCJ，Boelens HFM，et al．Analysis of longitudinalmetabolomics data．Bioinformatics，2004，20（15）：2438-2446．

［13］Jansen JJ，van Dam NM，Hoefsloot HCJ，etal．Crossfitanalysis：a novelmethod to characterize the dynamics of induced plant responses．BMC Bioinformatics，2009，10（1）：1．

［14］Ljung L．System identification．New Yersey：Prentice Hall，1987，163-173．

［15］Qin SJ，McAvoy TJ．Nonlinear FIR modeling via a neural net PLS approach．Comput Chem Eng，1996，20（2）：147-159．

［16］Dong Y，Qin SJ．Dynamic-Inner Partial Least Squares for Dynamic Data Modeling．IFAC-PapersOnLine，2015，48（8）：117-122．

［17］Kaspar MH，Ray WH．Dynamic PLSmodelling for process control．Chemical Engineering Science，1993，48（20）：3447-3461．

［18］Wold S，Kettaneh N，Fridén H，et al．Modelling and diagnostics of batch processesand analogous kinetic experiments．Chemometricsand intelligent laboratory systems，1998，44（1）：331-340．

［19］Antti H，Bollard ME，Ebbels T，et al．Batch statistical processing of 1H NMR-derived urinary spectral data．JChemom，2002，16（8-10）：461-468．

［20］Rantalainen M，Cloarec O，Ebbels TMD，et al．Piecewisemultivariate modelling of sequentialmetabolic profiling data．BMC Bioinformatics，2008，9（1）：105．

［21］Searle SR．Linear Models．New York：JW iley＆Sons，1971．

［22］Storey JD，Xiao W，Leek JT，et al．Significance analysis of time coursem icroarray experiments．Proc Natl Acad Sci USA，2005，102（36）：12837-12842．

［23］Conesa A，Nueda M J，Ferrer A，et al．maSigPro：amethod to identify significantly differential expression profiles in time-coursem icroarray experiments．Bioinformatics，2006，22（9）：1096-1102．

［24］Mardia KV，Kent JT，Bibby JM．Multivariate analysis．New York：Academ ic Press，1980．

［25］Sm ilde AK，Jansen JJ，Hoefsloot HCJ，et al．ANOVA-simultaneous component analysis（ASCA）：a new tool for analyzing designed metabolom ics data．Bioinformatics，2005，21（13）：3043-3048．

［26］Timmerman ME，Kiers HAL．Four simultaneous componentmodels for the analysis ofmultivariate time series from more than one subject to model intraindividual and interindividual differences．Psychometrika，2003，68（1）：105-121．

［27］Harrington PdB，Vieira NE，Espinoza J，et al．Analysis of varianceprincipal component analysis：A soft tool for proteom ic discovery．A-nal Chim Acta，2005，544（1）：118-127．

［28］de Noord OE，Theobald EH．Multilevel component analysis and multilevel PLS of chem ical process data．J Chemom，2005，19（5-7）：301-307．

［29］Jansen JJ，Bro R，Hoefsloot HCJ，et al．PARAFASCA：ASCA combined w ith PARAFAC for the analysis ofmetabolic fingerprinting data．JChemom，2008，22（2）：114-121．

［30］Sm ilde AK，Hoefsloot H，Westerhuis J．The geometry of ASCA．J Chemom，2008，22（8）：464-471．

［31］Mei Y，Kim SB，Tsui K．Linear-mixed effectsmodels for feature selection in high-dimensionalNMR spectra．Expert Systemswith Applications，2009，36（3）：4703-4708．

［32］Tipping ME，Bishop CM．Probabilistic principal component analysis．Journal of the Royal Statistical Society：Series B（Statistical Methodology），1999，61（3）：611-622．

［33］Jacquier E，Polson NG，Rossi PE．Bayesian analysis of stochastic volatility models．Journal of Business＆Economic Statistics，2002，20（1）：69-87．

［34］Tai YC，Speed TP．A multivariate empirical Bayes statistic for replicated microarray time course data．The Annals of Statistics，2006，34（5）：2387-2412．

［35］Xia J，Mandal R，Sinelnikov IV，et al．MetaboAnalyst2．0-a comprehensive server for metabolomic data analysis．Nucleic Acids Res，2012，40（W1）：W127-W133．

［36］Berk M，Ebbels T，Montana G．A statistical framework for biomarker discovery in metabolomic time course data．Bioinformatics，2011，27（14）：1979-1985．

［37］Wang J，Liu P，She MFH，et al．Biomedical time series clustering based on non-negative sparse coding and probabilistic topic model．ComputMethods Programs Biomed，2013，111（3）：629-641．

［38］Wang K，Ng SK，M cLachlan GJ．Clustering of time-course gene expression profiles using normal mixture models with autoregressive random effects．BMC Bioinformatics，2012，13（1）：1．

（責任編輯：郭海強）

*國家自然基金（81502889）；黑龍江省自然基金重點項目（ZD201314）

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