游懿君，梁丹丹，陳天璐

系統(tǒng)生物學是在細胞、組織、器官和生物體整體水平上研究結(jié)構(gòu)和功能各異的生物分子及其相互作用，并通過計算生物學來定量闡明和預(yù)測生物功能、表型和行為的科學[1]?；?、蛋白、轉(zhuǎn)錄、代謝及元基因等多組學整合研究是系統(tǒng)生物學中常用的研究手段?？缃M學研究不僅可了解機體的基因構(gòu)成，而且對于下游表達產(chǎn)物如蛋白、代謝物也有了一個充分的了解[1-2]。代謝組學是系統(tǒng)生物學領(lǐng)域一門新興的學科，是利用各儀器平臺(氣相/液相-色譜聯(lián)用和核磁共振等)和大數(shù)據(jù)信息挖掘策略觀察生物體系中的大量代謝產(chǎn)物在不同狀態(tài)下變化規(guī)律的科學[3]。元基因組學通過直接從環(huán)境樣品中提取全部微生物的DNA,構(gòu)建宏基因組文庫，利用基因組學的研究策略研究環(huán)境樣本所包含的全部微生物的遺傳組成及其群落功能。越來越多的研究表明，人體內(nèi)的代謝物水平波動與腸道內(nèi)菌群的變化間存在著一定的關(guān)聯(lián)[4-5]。由于人體是一個系統(tǒng)整體，系統(tǒng)性的研究顯得尤為必要[2]，而代謝組學和元基因組學2大組學平臺滿足了整體系統(tǒng)性分析的需求，越來越多的應(yīng)用到了菌群和代謝物的相關(guān)研究中。

相關(guān)性分析是指對2個或多個具備相關(guān)性的變量元素進行分析，從而衡量2個變量因素的相關(guān)密切程度[6]。利用相關(guān)性分析進行關(guān)聯(lián)變量的篩選是一種常用的科學研究策略。Pearson積相關(guān)[7]和Spearman秩相關(guān)[8]是2種經(jīng)典的相關(guān)性分析方法。隨著相關(guān)分析的廣泛應(yīng)用，在經(jīng)典的相關(guān)性方法基礎(chǔ)上發(fā)展延伸出了很多的適用于不同數(shù)據(jù)類型的方法。比如，SparCC[9]和CCLasso[10]是針對組分數(shù)據(jù)內(nèi)部變量的相關(guān)性研究而產(chǎn)生的方法，能夠更加精準的在該類數(shù)據(jù)集中找到關(guān)聯(lián)對。元基因組數(shù)據(jù)屬于組分數(shù)據(jù)的一種[11-12]，反映的是物質(zhì)的相對豐度信息而非絕對豐度信息，具有稀疏性和內(nèi)部關(guān)聯(lián)度高等特征。Friedman[9]和Fang等[10]指出，采用經(jīng)典相關(guān)分析(Pearson，Spearman)方法對其進行研究常常產(chǎn)生不合理甚至錯誤的結(jié)果。此外，尚未見有報道嘗試采用SparCC和CCLasso對非組分數(shù)據(jù)進行分析。代謝組學數(shù)據(jù)和元基因組學數(shù)據(jù)的數(shù)據(jù)特點和結(jié)構(gòu)也不盡相同。代謝組學數(shù)據(jù)是基于“質(zhì)譜”生成的信息，其變量是連續(xù)的且與具體的代謝物的濃度成比例；元基因組學數(shù)據(jù)是一種基于“計數(shù)”的信息，變量是離散的且與基因片段存在的數(shù)量成比例。2種不同特性組學數(shù)據(jù)的相關(guān)分析對方法的選擇提出了新的要求。

我們以代謝組學和元基因組學間的關(guān)聯(lián)研究為目標，選用4種典型的相關(guān)分析方法(Pearson，Spea-rman，SparCC和CCLasso)，采用仿真數(shù)據(jù)集和真實數(shù)據(jù)集，對這4種相關(guān)分析方法進行檢測和對比，研究了他們的相似性和獨特性，為跨組學數(shù)據(jù)相關(guān)分析方面提供方法學支持。

1 方法

1.1 仿真數(shù)據(jù)集建立 利用R軟件“MASS”包中的“mvrnorm”函數(shù)，生成2個含有21個變量(V)，200個樣本量(S)的數(shù)據(jù)矩陣A(a1-a21)和B(b1-b21)。其中代表代謝物的矩陣A服從正態(tài)分布，代表菌群的矩陣B服從對數(shù)正態(tài)分布。指定矩陣A與矩陣B中的21個相關(guān)對(a1-b1,a2-b2,a3-b3，…，a21-b21)的相關(guān)系數(shù)分別為-1,-0.9,-0.8，…,0.8,0.9,1。統(tǒng)計學上通常認為相關(guān)系數(shù)小于0.3是不存在相關(guān)或相關(guān)程度較弱，易產(chǎn)生不符合真實情況的結(jié)果，所以我們將相關(guān)系數(shù)大于或等于0.3的16個相關(guān)對作為有效相關(guān)對(編號第1～8和第14～21的相關(guān)對)。

1.2 腦代謝組研究數(shù)據(jù) 由于“菌群-腸-腦軸”的研究逐漸成為新的研究熱點，諸多研究證實腦中的代謝物與腸道菌群間存在著密切的相關(guān)性[13-15]。我們以42只Wistar大鼠腦組織代謝物的代謝組學數(shù)據(jù)和大鼠小腸內(nèi)容物中腸道菌群的元基因組學數(shù)據(jù)(16SrRNA)為例進行研究。所有的動物實驗都在上海交通大學實驗動物中心進行。原始的代謝組學數(shù)據(jù)包含了359個代謝物，原始的元基因組學數(shù)據(jù)包含了3 421個操作分類單元。我們對原始數(shù)據(jù)進行進一步的優(yōu)化：將代謝組學數(shù)據(jù)按照代謝物的種類進行合并，元基因組學數(shù)據(jù)操作分類單元表按分類學“門”水平進行合并。最終，得到30種代謝物和18個門的數(shù)據(jù)。

1.3 相關(guān)分析方法

1.3.1 Pearson積相關(guān)與Spearman秩相關(guān) Pearson積相關(guān)是Pearson[7]在1896年創(chuàng)立的一種經(jīng)典相關(guān)分析方法(公式1)。在此基礎(chǔ)上后續(xù)延伸出了Spearman秩相關(guān)[8]。Spearman從公式上來看是一種非參的Pearson相關(guān)，在計算相關(guān)系數(shù)時，不體現(xiàn)計算的具體數(shù)值，而是將具體數(shù)值按秩次進行排序，利用秩次進行相關(guān)系數(shù)的計算。2種方法皆在R中調(diào)用“cor”函數(shù)實現(xiàn)，選擇method=“pearson”或“spearman”。

(1)

1.3.2 SparCC與CCLasso SparCC和CCLasso是主要用于分析組分數(shù)據(jù)間相關(guān)關(guān)系的一種非參的相關(guān)分析方法，是前期經(jīng)典相關(guān)分析方法的延伸。SparCC方法的大體思路是：將組分數(shù)據(jù)進行Aitchison對數(shù)轉(zhuǎn)換，計算轉(zhuǎn)換后的線性Pearson相關(guān)并得到相關(guān)系數(shù)[9]。CCLasso則是在對數(shù)轉(zhuǎn)換后對組分進行基于罰函數(shù)的最小二乘法(即Lasso)從而得到相關(guān)系數(shù)[10]。SparCC和CCLasso方法均可在R軟件中實現(xiàn)，方法代碼可從https://github.com/huayingfang/CCLasso中獲得[10]。

抗胃癌植物類中藥藥味以苦、甘、辛為主，豐度分別為0.612、0.518、0.271；主要涉及品種包括白花蛇舌草、丹參、纈草、柴胡等，詳見表3。臨床在選擇組方入藥時應(yīng)結(jié)合患者臨床癥狀及個體特征充分考慮上述藥味特點。

本研究的所有相關(guān)方法的實現(xiàn)皆在R軟件中進行。

2 結(jié)果與討論

將4種相關(guān)方法對仿真數(shù)據(jù)集和腦代謝組真實數(shù)據(jù)集進行相關(guān)分析，并將結(jié)果進行統(tǒng)計分析和可視化處理，用以從不同角度衡量對比各方法的性能。

從圖1A中可以得出，16個指定相關(guān)對經(jīng)4種相關(guān)方法計算得到的r值(-1～-0.3，0.3～1)大體一致。圖1B中可以得出CCLasso的誤差百分比大于其他3種方法(22.95%)。Pearson和Spearman方法得出的誤差百分比較小(Pearson:6.23%;Spearman:5.14%)。圖1的結(jié)果說明了CCLasso方法存在較大的誤差，計算出的相關(guān)系數(shù)與真實值差距較大。這可能是由于CCLasso算法復(fù)雜度高，對數(shù)據(jù)結(jié)構(gòu)要求較嚴格。另一個原因可能是由于前期仿真數(shù)據(jù)集的構(gòu)建是根據(jù)Pearson系數(shù)來設(shè)計的(協(xié)方差矩陣公式是Pearson相關(guān)的一種變換)，所以Pearson和Spearman相對來說得到的結(jié)果誤差百分比偏小，結(jié)果偏好。盡管如此，相比之下CCLasso的準確度比SparCC稍差。

A:不同方法對指定r值的16個相關(guān)對進行分析得到的相關(guān)系數(shù)(r) 所有P值均小于0.05

B:4種相關(guān)方法所計算出的相關(guān)系數(shù)的誤差百分比(均值±標準誤)。 誤差百分比=100%×|計算值-設(shè)定值|/設(shè)定值。4組間ANOVA的P<0.001， *表示兩兩比較時(Tukey’s方法)，該方法與CCLasso比較的P<0.05圖1 4種相關(guān)方法基于仿真數(shù)據(jù)集中指定了相關(guān)系數(shù) 的16個變量對的計算結(jié)果

進一步采用4種相關(guān)方法對矩陣A、B中未指定相關(guān)系數(shù)的變量對(變量如a1與b2,b3,…,b21等)進行分析，其相關(guān)結(jié)果如圖2所示。由r值聚類熱圖(聚類分析中使用歐幾里得距離計算變量之間的距離)(圖2A)可得,Pearson與Spearman聚為一類，SparCC與CCLasso聚為一類。該結(jié)果表明方法Pearson與Spearman得到的結(jié)果更為相似，SparCC與CCLasso的結(jié)果更為相似。而從數(shù)學公式上比較，這2組方法計算相關(guān)系數(shù)的公式也更為相近。圖2B顯示的是各種方法找到的顯著相關(guān)對的數(shù)目。我們發(fā)現(xiàn)，CCLasso找到的相關(guān)對(P<0.05)數(shù)目最少(2)，SparCC最多(23)，Spearman(21)和Pearson(18)居中。這些未指定相關(guān)系數(shù)的變量對可相當于隨機變量，理論上這些變量間不存在相關(guān)性。該結(jié)果說明CCLasso最為保守，經(jīng)典方法的性能居中，而采用SparCC的假陽性風險相對較高。在必要時，可采用各種方法對P值進行適當校正或采用更嚴格的閾值，以降低假陽性。

A:基于相關(guān)系數(shù)r值的層次聚類結(jié)果

B:顯著相關(guān)(P<0.05)的個數(shù)圖2 4種相關(guān)方法對仿真數(shù)據(jù)集中未指定相關(guān)系數(shù)的 變量對(21×20=420個)的計算結(jié)果

最后，采用真實數(shù)據(jù)對各種方法的性能進行進一步比較和評估。根據(jù)圖3A所示，對于同樣的數(shù)據(jù)集，CCLasso計算得到的相關(guān)系數(shù)r值相對最小，SparCC計算所得的r值最高，Spearman和Pearson介于兩者之間。圖3C顯示的是顯著(P<0.05)相關(guān)對數(shù)目。與仿真數(shù)據(jù)集和圖3A結(jié)果一致的是，CCLasso找到的相關(guān)對數(shù)目最少，SparCC最多，Spearman和Pearson結(jié)果介于兩者之間。由此我們判斷CCLasso方法相對保守，對數(shù)據(jù)的要求比較嚴格。此外，該方法的計算耗時也較其他方法長。這可能是因為其內(nèi)部算法比較復(fù)雜引起的[10]。再次采用相關(guān)系數(shù)r進行的層次聚類分析(圖3B)表明，方法Pearson與Spearman得到的結(jié)果更為相似，方法SparCC與CCLasso的結(jié)果更為相似。該結(jié)果與仿真數(shù)據(jù)集聚類結(jié)果一致。

從仿真數(shù)據(jù)集和真實數(shù)據(jù)集的實驗結(jié)果來看，我們認為傳統(tǒng)的相關(guān)分析方法在跨組學之間的相關(guān)分析中更具有優(yōu)勢。而Spearman由于適合非線性相關(guān)，建議其為首選方法。

A:相關(guān)系數(shù)r值的范圍。橫坐標為4種相關(guān)方法， 縱坐標為相關(guān)系數(shù)r值(均值±標準差)。 4組間ANOVA的P<0.001。兩兩比較(Tukey’s)的所有P<0.05

B:基于相關(guān)系數(shù)r值的層次聚類和熱圖

C:顯著(P<0.05)相關(guān)對數(shù)目圖3 4種相關(guān)方法基于真實數(shù)據(jù)集的計算結(jié)果

3 總結(jié)與展望

基于以上仿真和真實數(shù)據(jù)集的結(jié)果，可以發(fā)現(xiàn)4種相關(guān)分析方法計算得到的r值和顯著相關(guān)對的數(shù)量都存在類似的規(guī)律：SparCC>Spearman和Pearson>CCLasso。CCLasso誤差百分比較大，而得到的顯著相關(guān)對的數(shù)目最少。相比較來說，該方法的結(jié)果更容易出現(xiàn)假陰性(即具有相關(guān)關(guān)系的2個變量關(guān)系，在CCLasso較為嚴格的算法下，結(jié)果往往是不具有相關(guān)性)。相反的，SparCC更傾向于給出假陽性結(jié)果(即不具有相關(guān)關(guān)系的2個變量關(guān)系，在SparCC較為寬松的算法下，結(jié)果往往是具有相關(guān)性)；經(jīng)典相關(guān)方法(Pearson和Spearman)的表現(xiàn)較為穩(wěn)定平和，Spearman結(jié)果稍優(yōu)于Pearson。部分原因可能是其更適合于分析非線性相關(guān)關(guān)系。本研究結(jié)論仍需要更多實驗數(shù)據(jù)的進一步驗證。

【參考文獻】

[1]楊勝利.系統(tǒng)生物學研究進展[J].中國科學院院刊,2004,19(1):31-34.

[2]Kitano H.Systems biology:a brief overview[J].Science,2002,295(5560):1662-1664.

[3]Bujak R,Struck-Lewicka W,Markuszewski MJ,et al.Metabolomics for laboratory diagnostics[J].J Pharm Biomed Anal,2015,113:108-120.

[4]Liu R,Hong J,Xu X,et al.Gut microbiome and serum metabolome alterations in obesity and after weight-loss intervention[J].Nat Med,2017,23(7):859-868.

[5]Zheng X,Xie G,Zhao A,et al.The footprints of gut microbial-mammalian co-metabolism[J].J Proteome Res,2011,10(12):5512-5522.

[6]Hauke J,Kossowski T.Comparison of values of Pearson’s and Spearman’s correlation coefficients on the same sets of data[J].Quaest Geo,2011,30(2):87-93.

[7]Pearson K.Mathematical contributions to the theory of evolution.III.Regression, heredity,and panmixia[J].Philos T Roy Soc A,1896,187:253-318.

[8]Spearman C.The proof and measurement of association between two things[J].Am J Psychol,1904,15(1):72-101.

[9]Friedman J,Alm EJ.Inferring correlation networks from genomic survey data[J].Plos Comput Biol,2012,8(9):e1002687.

[10]Fang H,Huang C,Zhao H,et al.CCLasso: correlation inference for compositional data through Lasso[J].Bioinformatics,2015,31(19):3172-3180.

[11]Tsilimigras MC,Fodor AA.Compositional data analysis of the microbiome:fundamentals,tools,and challenges[J].Ann Epidemiol,2016,26(5):330-335.

[12]Weiss S,Van Treuren W,Lozupone C,et al.Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[J].ISME J,2016,10(7):1669-1681.

[13]Zheng X,Chen T,Zhao A,et al.The brain metabolome of male rats across the lifespan[J].Sci Rep,2016,6:24125.

[14]Mayer EA.Gut feelings:the emerging biology of gut-brain communication[J].Nat Rev Neurosci,2011,12(8):453-466.

[15]Ridaura V,Belkaid Y.Gut microbiota:the link to your second brain[J].Cell,2015,161(2):193-194.