基于gCoda方法的HMP成分?jǐn)?shù)據(jù)分析

徐平峰，牛彩虹，孫 萌，王樹(shù)達(dá)

(長(zhǎng)春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，吉林 長(zhǎng)春 130012)

0 引 言

自然環(huán)境中微生物隨處可見(jiàn)，在人體健康中扮演著重要角色。分析微生物群落間的關(guān)系可以探索微生物對(duì)我們生活環(huán)境的影響，而研究微生物的一個(gè)重要目標(biāo)是推斷微生物間的相互作用網(wǎng)絡(luò)，這就涉及到了圖模型。孫聚波等[1]在前人的基礎(chǔ)上對(duì)圖模型的基本概念進(jìn)行了簡(jiǎn)要回顧?，F(xiàn)有研究針對(duì)微生物的相互作用網(wǎng)絡(luò)提出了幾種方法;Friedman等[2]基于成分?jǐn)?shù)據(jù)的對(duì)數(shù)比引入潛變量的概念，提出了Sparcc(Sparse Correlations for Compositional data)的近似方法來(lái)推斷稀疏假設(shè)下的相關(guān)矩陣，但其沒(méi)有考慮到成分?jǐn)?shù)據(jù)中誤差的影響，估計(jì)出的協(xié)方差矩陣也不能保證其正定性；Fang等[3]提出了CCLasso(Correlation inference for Compositional data through Lasso)，這種方法將加權(quán)最小二乘和l1懲罰結(jié)合起來(lái)推斷相關(guān)矩陣，與Sparcc相比，其性能更好，但這些方法都沒(méi)有統(tǒng)計(jì)假設(shè)可以保證能有效進(jìn)行;Huaying Fang等[4]提出了一種新的懲罰最大似然方法gCoda(generation mechanism of compositional data)，從觀察到的成分?jǐn)?shù)據(jù)logistic正態(tài)分布推斷微生物之間的稀疏直接相互作用網(wǎng)絡(luò)，估計(jì)潛變量的逆協(xié)方差的稀疏結(jié)構(gòu)。

文中采用gCoda方法分析HMP(NIH Human Microbiome Project)的數(shù)據(jù)，刻畫(huà)人體中由基因組數(shù)據(jù)產(chǎn)生的成分?jǐn)?shù)據(jù)的相關(guān)性，進(jìn)而分析微生物相互間的作用關(guān)系，為研究微生物對(duì)人類(lèi)健康和疾病等方面提供一些依據(jù)。

1 成分?jǐn)?shù)據(jù)

Aitchison[5]首次提出對(duì)成分?jǐn)?shù)據(jù)進(jìn)行l(wèi)og-ratio變換，使得對(duì)成分?jǐn)?shù)據(jù)的研究更進(jìn)一步。假設(shè)y=(y1,y2,…,yp)T是p維隨機(jī)計(jì)數(shù)潛變量，x=(x1,x2,…,xp)T能夠被直接觀測(cè)，且x與y的關(guān)系為

(1)

j=1,2,…,p,

lnx=lny-1plnw,

(2)

式中：1p----元素全為1的p維列向量。

2 gCoda

Fang H等[3]提出的gCoda方法中，假設(shè)成分?jǐn)?shù)據(jù)服從logistic正態(tài)分布，可以將我們的推斷轉(zhuǎn)換為正態(tài)分布的逆協(xié)方差問(wèn)題；假設(shè)相互作用網(wǎng)絡(luò)是稀疏的，這可以解決由成分性引起的不確定性問(wèn)題。

假設(shè)lny～Np(μ,Ω-1)，則(lnw,x)的聯(lián)合密度分布為

(3)

式中：M=lnx+1plnw-μ。

令變換矩陣

式中：Ep----p維單位矩陣。

對(duì)f(lnw,x)積分,得到f(x)的分布

式中:N=F0lnx-F0μ。

樣本X=(x1,x2,…,xp)是獨(dú)立且同分布的，其負(fù)對(duì)數(shù)似然L(μ,Ω)為

(4)

其中

假設(shè)E(F0lnx)的估計(jì)為E(F0μ)，得到Ω的負(fù)對(duì)數(shù)似然函數(shù)為

(5)

其中

處理稀疏協(xié)方差最常用的方法是給損失函數(shù)加一個(gè)懲罰項(xiàng)，所以用下面的函數(shù)作為我們的目標(biāo)函數(shù)

f(Ω)=L(Ω)+λn‖Ω‖1,

(6)

其中

式中：λn----調(diào)節(jié)參數(shù)，λn>0。

gCoda的目的是找到懲罰似然的最大似然估計(jì)，即

(7)

由于式(7)是非凸的，所以用MM算法來(lái)迭代(Lange等[6])求解目標(biāo)函數(shù)的最小值。當(dāng)算法進(jìn)行到第k步時(shí)，f(Ω)的替代函數(shù)為g(Ω)，即

gk(Ωk)=-ln|Ω|+tr(Ω(Ep-O)S0(Ep-P))+

λn‖Ω‖1,

得出

gk(Ωk)=f(Ωk)。

gk(Ωk)可以被看做是一個(gè)標(biāo)準(zhǔn)的glasso問(wèn)題，因?yàn)?/p>

tr(ΩSk)+λn‖Ω‖1,

(8)

其中

最終把MM算法的優(yōu)化問(wèn)題轉(zhuǎn)化成通過(guò)坐標(biāo)軸下降法可以解決的glasso[7]問(wèn)題，同時(shí)用BIC(Bayesian Information Criterion)來(lái)選擇懲罰參數(shù)。

3 數(shù)值模擬及實(shí)例分析

3.1 數(shù)值模擬

使用成分?jǐn)?shù)據(jù)而非計(jì)數(shù)數(shù)據(jù)做如下模擬。考慮服從logistic正態(tài)分布的成分?jǐn)?shù)據(jù)

X=(x1,x2,…,xp)，

其中:

xj=(x1j,x2j,…,xnj)T，

i=1,2,…,n，j=1,2,…,p，

lny～Np(μ,Ω-1)。

lny的均值向量μ從均勻分布[-0.5,0.5]p中產(chǎn)生,同時(shí),用無(wú)標(biāo)度圖作為lny的協(xié)方差逆陣。

無(wú)標(biāo)度圖(Scale-free)使用B-A算法[8]建立一個(gè)無(wú)標(biāo)度圖。從單個(gè)頂點(diǎn)開(kāi)始，然后一個(gè)接一個(gè)地添加p-1個(gè)頂點(diǎn)。每次新增一個(gè)節(jié)點(diǎn)時(shí)，新節(jié)點(diǎn)都與三個(gè)隨機(jī)選擇的舊節(jié)點(diǎn)相連，這些舊節(jié)點(diǎn)被選擇的概率與當(dāng)前圖中每個(gè)節(jié)點(diǎn)的自由度有關(guān)。邊緣強(qiáng)度在區(qū)間[-0.8,-0.6]和[0.6,0.8]的均勻分布中隨機(jī)產(chǎn)生。為確保Ω是一個(gè)正定矩陣，我們把Ω的對(duì)角線元素都設(shè)置為5。

為評(píng)估gCoda的性能，將tpr(真陽(yáng)性比率)、ppv(陽(yáng)性預(yù)測(cè)率)和mcc(馬修斯相關(guān)系數(shù))作為評(píng)價(jià)指標(biāo)。其中:

式中：tp，tn，fp，fn----分別為真陽(yáng)性、真陰性、假陽(yáng)性、假陰性的數(shù)量。

mcc的取值范圍為[-1,1]，當(dāng)?shù)貌坏絤cc的值時(shí)，將mcc的值記為-1。

設(shè)置變量個(gè)數(shù)p=50，樣本量n=(100,200,500)，且對(duì)每個(gè)模型都進(jìn)行50次模擬，計(jì)算tpr、ppv和mcc的值來(lái)評(píng)估gCoda的估計(jì)性能。

在不同樣本量下，gCoda方法對(duì)無(wú)標(biāo)度圖進(jìn)行估計(jì)，得到三個(gè)評(píng)價(jià)指標(biāo)的均值與方差見(jiàn)表1。

表1 不同樣本量下，gCoda的性能(標(biāo)準(zhǔn)差)

由表1可以看出,隨著樣本量的增加，tpr和mcc的均值都在增加，即樣本量越大,越可以估計(jì)得到更多的真邊，估計(jì)所得到的圖與真實(shí)圖也越接近。

gCoda估計(jì)邊數(shù)量的直方圖如圖1所示。

圖1 gCoda估計(jì)邊數(shù)量的直方圖

圖中，黑色代表tp的均值，灰色代表fp的均值，tp+fp是估計(jì)圖中邊的數(shù)量，而圖中的橫線表示原始圖中存在邊的數(shù)量。每幅圖中橫軸表示樣本量的大?。豢v軸代表tp+fp的數(shù)量。以上結(jié)果基于50次模擬得出的均值。

通過(guò)圖1可以看到，隨著樣本量的增加，tp和fp的數(shù)量均增大，且當(dāng)樣本量為500時(shí)，tp的個(gè)數(shù)已經(jīng)非常接近真實(shí)圖中邊的個(gè)數(shù)。

3.2 實(shí)證分析

OTUs是生物信息分析中一種常見(jiàn)的計(jì)數(shù)數(shù)據(jù)，文中所用OUT計(jì)數(shù)是不同人體的頰粘膜、喉嚨等18個(gè)身體部位的樣本。原始數(shù)據(jù)以及樣本說(shuō)明可以從網(wǎng)址http://hmpdacc.org/HMQCP/上獲取。文中刪除數(shù)據(jù)中OTUs為0的個(gè)數(shù)>60%的樣本。值得注意的是，式(2)對(duì)成分?jǐn)?shù)據(jù)進(jìn)行對(duì)數(shù)運(yùn)算，所以，文中對(duì)計(jì)數(shù)數(shù)據(jù)加了一個(gè)值為0.5的偽計(jì)數(shù)作為分析所用數(shù)據(jù)。除此之外，文中將分類(lèi)不明確的OTUs剔除，并且按照每行、每列OTUs總計(jì)數(shù)不小于500篩選一部分樣本，最終保留了2 744個(gè)樣本，并選擇其中60個(gè)變量進(jìn)行數(shù)據(jù)分析。

用BIC準(zhǔn)則選擇得到的最優(yōu)模型，如圖2所示。

圖2 gCoda估計(jì)圖

由圖2可以看到，編號(hào)54、10、22與其他所有變量都沒(méi)有關(guān)系，而位于圖中左下角的變量間具有復(fù)雜的相關(guān)性。同屬細(xì)菌之間不僅具有相互關(guān)系，不同屬細(xì)菌之間也可能存在關(guān)系，這可能與其生活環(huán)境有關(guān)。如編號(hào)11、12和44都是Bacteroides屬，而編號(hào)34屬于Ruminococcus屬，編號(hào)48是Prevotella屬，編號(hào)44與不同屬的34、48具有相互關(guān)系，而與同屬的11、12有關(guān)系。

4 結(jié) 語(yǔ)

用gCoda方法對(duì)成分?jǐn)?shù)據(jù)進(jìn)行數(shù)值模擬，并用tpr、ppv和mcc值作為評(píng)估此種算法的性能指標(biāo)，同時(shí)對(duì)NIH Human Microbiome Project數(shù)據(jù)進(jìn)行分析，得到微生物間的相互作用網(wǎng)絡(luò)，可能對(duì)人體微生物與人體健康的關(guān)系有一定幫助。