基于gCoda方法的HMP成分數據分析

徐平峰，牛彩虹，孫 萌，王樹達

(長春工業大學 數學與統計學院，吉林 長春 130012)

0 引 言

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

文中采用gCoda方法分析HMP(NIH Human Microbiome Project)的數據，刻畫人體中由基因組數據產生的成分數據的相關性，進而分析微生物相互間的作用關系，為研究微生物對人類健康和疾病等方面提供一些依據。

1 成分數據

Aitchison[5]首次提出對成分數據進行log-ratio變換，使得對成分數據的研究更進一步。假設y=(y1,y2,…,yp)T是p維隨機計數潛變量，x=(x1,x2,…,xp)T能夠被直接觀測，且x與y的關系為

(1)

j=1,2,…,p,

lnx=lny-1plnw,

(2)

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

2 gCoda

Fang H等[3]提出的gCoda方法中，假設成分數據服從logistic正態分布，可以將我們的推斷轉換為正態分布的逆協方差問題；假設相互作用網絡是稀疏的，這可以解決由成分性引起的不確定性問題。

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

(3)

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

令變換矩陣

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

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

式中:N=F0lnx-F0μ。

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

(4)

其中

假設E(F0lnx)的估計為E(F0μ)，得到Ω的負對數似然函數為

(5)

其中

處理稀疏協方差最常用的方法是給損失函數加一個懲罰項，所以用下面的函數作為我們的目標函數

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

(6)

其中

式中：λn----調節參數，λn>0。

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

(7)

由于式(7)是非凸的，所以用MM算法來迭代(Lange等[6])求解目標函數的最小值。當算法進行到第k步時，f(Ω)的替代函數為g(Ω)，即

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

λn‖Ω‖1,

得出

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

gk(Ωk)可以被看做是一個標準的glasso問題，因為

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

(8)

其中

最終把MM算法的優化問題轉化成通過坐標軸下降法可以解決的glasso[7]問題，同時用BIC(Bayesian Information Criterion)來選擇懲罰參數。

3 數值模擬及實例分析

3.1 數值模擬

使用成分數據而非計數數據做如下模擬。考慮服從logistic正態分布的成分數據

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中產生,同時,用無標度圖作為lny的協方差逆陣。

無標度圖(Scale-free)使用B-A算法[8]建立一個無標度圖。從單個頂點開始，然后一個接一個地添加p-1個頂點。每次新增一個節點時，新節點都與三個隨機選擇的舊節點相連，這些舊節點被選擇的概率與當前圖中每個節點的自由度有關。邊緣強度在區間[-0.8,-0.6]和[0.6,0.8]的均勻分布中隨機產生。為確保Ω是一個正定矩陣，我們把Ω的對角線元素都設置為5。

為評估gCoda的性能，將tpr(真陽性比率)、ppv(陽性預測率)和mcc(馬修斯相關系數)作為評價指標。其中:

式中：tp，tn，fp，fn----分別為真陽性、真陰性、假陽性、假陰性的數量。

mcc的取值范圍為[-1,1]，當得不到mcc的值時，將mcc的值記為-1。

設置變量個數p=50，樣本量n=(100,200,500)，且對每個模型都進行50次模擬，計算tpr、ppv和mcc的值來評估gCoda的估計性能。

在不同樣本量下，gCoda方法對無標度圖進行估計，得到三個評價指標的均值與方差見表1。

表1 不同樣本量下，gCoda的性能(標準差)

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

gCoda估計邊數量的直方圖如圖1所示。

圖1 gCoda估計邊數量的直方圖

圖中，黑色代表tp的均值，灰色代表fp的均值，tp+fp是估計圖中邊的數量，而圖中的橫線表示原始圖中存在邊的數量。每幅圖中橫軸表示樣本量的大小；縱軸代表tp+fp的數量。以上結果基于50次模擬得出的均值。

通過圖1可以看到，隨著樣本量的增加，tp和fp的數量均增大，且當樣本量為500時，tp的個數已經非常接近真實圖中邊的個數。

3.2 實證分析

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

用BIC準則選擇得到的最優模型，如圖2所示。

圖2 gCoda估計圖

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

4 結 語

用gCoda方法對成分數據進行數值模擬，并用tpr、ppv和mcc值作為評估此種算法的性能指標，同時對NIH Human Microbiome Project數據進行分析，得到微生物間的相互作用網絡，可能對人體微生物與人體健康的關系有一定幫助。