動態(tài)代謝組學數(shù)據(jù)分析方法介紹*

哈爾濱醫(yī)科大學公共衛(wèi)生學院衛(wèi)生統(tǒng)計學教研室（150081） 王文佶 張秋菊 曲思楊 謝彪 高兵 孫琳 劉美娜

動態(tài)代謝組學數(shù)據(jù)分析方法介紹*

哈爾濱醫(yī)科大學公共衛(wèi)生學院衛(wèi)生統(tǒng)計學教研室（150081） 王文佶 張秋菊 曲思楊 謝彪 高兵 孫琳 劉美娜△

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

常用代謝組學數(shù)據(jù)分析方法

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

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

動態(tài)代謝組學數(shù)據(jù)分析方法

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

1．動態(tài)主成分分析

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

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

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

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

2．動態(tài)偏最小二乘法

這里介紹三種動態(tài)偏最小二乘法：

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

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

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

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

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

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

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

4．動態(tài)概率主成分分析（dynamic probabilistic principal components analysis，DPPCA）

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

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

DPPCA模型的構(gòu)建原理為：假設(shè)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模型。

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

5．其他方法：

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

小 結(jié)

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

在選擇一種動態(tài)分析模型時應(yīng)考慮實驗設(shè)計所涉及的時間點的個數(shù)、時間的間隔、所測得的代謝物的數(shù)量?？梢愿鶕?jù)研究設(shè)計和目的不同選用加權(quán)PCA，局部PCA批次模型從不同角度解釋時間帶來的變異性。若想檢驗數(shù)據(jù)不同來源的變異隨時間的變化趨勢可選用ASCA模型；若研究關(guān)注隨時間改變不同處理組的差異代謝產(chǎn)物可選用DPPCA和混合線性模型；另外，多時點的動態(tài)代謝組學數(shù)據(jù)慎重選用動態(tài)PCA避免因時間點增多丟失的信息逐漸增加。

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

［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，F(xiàn)ridé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，F(xiàn)errer 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