4種相關(guān)分析方法在菌群和代謝物相關(guān)研究中的初步比較

游懿君，梁丹丹，陳天璐

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

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

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

1 方法

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

1.2 腦代謝組研究數(shù)據(jù) 由于“菌群-腸-腦軸”的研究逐漸成為新的研究熱點(diǎn)，諸多研究證實(shí)腦中的代謝物與腸道菌群間存在著密切的相關(guān)性[13-15]。我們以42只Wistar大鼠腦組織代謝物的代謝組學(xué)數(shù)據(jù)和大鼠小腸內(nèi)容物中腸道菌群的元基因組學(xué)數(shù)據(jù)(16SrRNA)為例進(jìn)行研究。所有的動(dòng)物實(shí)驗(yàn)都在上海交通大學(xué)實(shí)驗(yàn)動(dòng)物中心進(jìn)行。原始的代謝組學(xué)數(shù)據(jù)包含了359個(gè)代謝物，原始的元基因組學(xué)數(shù)據(jù)包含了3 421個(gè)操作分類(lèi)單元。我們對(duì)原始數(shù)據(jù)進(jìn)行進(jìn)一步的優(yōu)化：將代謝組學(xué)數(shù)據(jù)按照代謝物的種類(lèi)進(jìn)行合并，元基因組學(xué)數(shù)據(jù)操作分類(lèi)單元表按分類(lèi)學(xué)“門(mén)”水平進(jìn)行合并。最終，得到30種代謝物和18個(gè)門(mén)的數(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從公式上來(lái)看是一種非參的Pearson相關(guān)，在計(jì)算相關(guān)系數(shù)時(shí)，不體現(xiàn)計(jì)算的具體數(shù)值，而是將具體數(shù)值按秩次進(jìn)行排序，利用秩次進(jìn)行相關(guān)系數(shù)的計(jì)算。2種方法皆在R中調(diào)用“cor”函數(shù)實(shí)現(xiàn)，選擇method=“pearson”或“spearman”。

(1)

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

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

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

2 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

3 總結(jié)與展望

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

【參考文獻(xiàn)】

[1]楊勝利.系統(tǒng)生物學(xué)研究進(jìn)展[J].中國(guó)科學(xué)院院刊,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.