一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法及裝置的制造方法
【專利摘要】本發(fā)明提供一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法及裝置，該方法包括：獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組測序樣本對應(yīng)的環(huán)境因素測量值；根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組測序樣本中各種微生物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型；采用最大后驗估計算法對分層貝葉斯模型進行學(xué)習(xí)，確定分層貝葉斯模型的目標(biāo)函數(shù)；對目標(biāo)函數(shù)進行優(yōu)化；采用優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)，從而提高在微生物關(guān)聯(lián)和微生物與環(huán)境因素關(guān)聯(lián)的預(yù)測任務(wù)中的準(zhǔn)確性和實用性。
【專利說明】
一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法及裝置
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及生物信息技術(shù)領(lǐng)域，尤其涉及一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法及裝置。
【背景技術(shù)】
[0002] 微生物之間的交互以及環(huán)境因素對微生物的影響是微生物學(xué)的關(guān)鍵的研究課題。 由于大部分的微生物都無法在實驗室培養(yǎng)，這給傳統(tǒng)的利用微生物培養(yǎng)來研究交互作用的 方法帶來了很大困難。宏基因組測序技術(shù)的發(fā)展為研究這些問題帶來了可能性，可以通過 采集微生物的樣本(提取16s rDNA)，然后進行高通量測序，通過樣本中微生物的數(shù)量變化 來推測微生物之間的關(guān)聯(lián)，以及環(huán)境因素對微生物的影響。隨著宏基因組測序技術(shù)的廣泛 應(yīng)用，各種環(huán)境的樣本變得可用，如土壤，海洋和湖泊，還有很多人的腸道。但是通過微生物 在樣本中的數(shù)量變化來推斷交互作用，仍然是一個很難的問題。通常的關(guān)聯(lián)推斷，指的是計 算微生物之間，微生物和環(huán)境之間的正負(fù)相關(guān)性，然后根據(jù)這些相關(guān)性再去推測真實的交 互作用，如互利共生、寄生和競爭，從而幫助理解微生物群落的動態(tài)。
[0003] 關(guān)聯(lián)推斷可以通過不同的統(tǒng)計量度量，只要這些統(tǒng)計量能夠顯示合理的關(guān)系。現(xiàn) 有的關(guān)聯(lián)推斷的方法可以分為兩大類。一種是計算兩兩之間的相關(guān)性，例如皮爾遜相關(guān)系 數(shù)和斯皮爾曼相關(guān)系數(shù)，可以計算兩個物種之間的相關(guān)性。Ruanet.al. (2006)年提出的 Local similarity analysis(LSA)計算的也是兩兩相關(guān)性，但是需要時間序列信息。第二 類是計算復(fù)雜相關(guān)性的方法，通過擬合一種微生物，和剩余微生物以及環(huán)境因素之間的關(guān) 系進行計算，大多是基于回歸的方法。第一種計算兩兩之間相關(guān)性的方法因為其簡單快速 的特點被生物學(xué)家廣泛采用。但是這種方法并不適合宏基因組測序數(shù)據(jù)，主要原因有兩點。 第一，計算兩兩相關(guān)性的方法并沒有得到真實的相關(guān)性，因為計算過程中存在組成成分偏 差，而這些方法忽略了數(shù)據(jù)的特點，采用的是適用于無約束數(shù)據(jù)的方法。具體來講，因為測 序得到的每個樣本的總reads數(shù)不同，一般會進行歸一化，得到OTU的相對豐度。在歸一化之 后，數(shù)據(jù)便不再獨立，例如變量剩余的變量&不再獨立，無論他們之間有無關(guān)聯(lián)，
[0004]
[0005] 這種組成成分偏差在環(huán)境中存在主導(dǎo)的微生物是會變得更加嚴(yán)重，這種現(xiàn)象其實 在海洋微生物中普遍存在。因此，考慮組成成分偏差的算法是需要的。另一方面，高通量測 序伴隨著一系列處理步驟，如樣本過濾、擴增和上機測序等，這些步驟都會導(dǎo)致得到的 reads數(shù)和樣本中原來的微生物數(shù)量出現(xiàn)偏差。這個特點也需要算法同時考慮。
[0006] 目前，出現(xiàn)了很多算法來解決組成成分偏差的問題，如CCREPE(Faust et.al.， 2012)，SparCC(Friedman and Alm，2012)，SPIEC-EASI(Kurtz et.al.,2015)和CCLasso (Fang et.al.，2015)。這些算法通過不同的思路處理組成成分偏差，但是均未考慮測序數(shù) 據(jù)本身的方差;另外，由于微生物會受到環(huán)境因素的影響，這些方法在估計微生物之間的關(guān) 聯(lián)時，并未考慮環(huán)境因素的調(diào)控。例如，如果兩個OTU有關(guān)聯(lián)如果是因為受到同一種環(huán)境因 素的調(diào)節(jié)，那么他們之間的關(guān)聯(lián)其實是間接的，如果不考慮環(huán)境因素的話，是無法做區(qū)分 的。

【發(fā)明內(nèi)容】

[0007] 鑒于上述問題，本發(fā)明提出了一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法及裝置，以解決現(xiàn)有 技術(shù)中因為組成成分偏差和測序數(shù)據(jù)自身的方差所導(dǎo)致的關(guān)聯(lián)推斷不準(zhǔn)確的問題。
[0008] 根據(jù)本發(fā)明的第一方面，提供了一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法，該方法包括：
[0009] 獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組測序樣本對應(yīng)的 環(huán)境因素測量值；
[0010] 根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組測序樣本中各種微生 物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型，所述分層貝葉斯模型用于對 數(shù)據(jù)產(chǎn)生過程、微生物之間的關(guān)聯(lián)以及微生物與環(huán)境因素測量值之間的關(guān)聯(lián)進行描述； [0011]采用最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分層貝葉斯模 型的目標(biāo)函數(shù)；
[0012] 對所述目標(biāo)函數(shù)進行優(yōu)化；
[0013] 采用目標(biāo)函數(shù)優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和微生物與環(huán)境 因素測量值之間的關(guān)聯(lián)。
[0014] 其中，所述根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組測序樣本 中各種微生物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型，包括：
[0015] 采用Dirichlet-Multinomial共輒分布模擬測序過程，對數(shù)據(jù)產(chǎn)生過程中的組成 成分偏差和數(shù)據(jù)方差進行第一建模；
[0016]采用Lognormal分布模擬微生物的豐度變化過程，對微生物之間的關(guān)聯(lián)以及微生 物與環(huán)境因素測量值之間的關(guān)聯(lián)進行第二建模；
[0017] 根據(jù)所述第一建模和第二建模的結(jié)果，建立基于LognormaI-D iri ch I et-Multinomial的分層貝葉斯模型。
[0018] 其中，所述采用最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分 層貝葉斯模型的目標(biāo)函數(shù)，包括：
[0019] 采用增加一范數(shù)懲罰項的最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確 定所述分層貝葉斯模型的目標(biāo)函數(shù)。
[0020] 其中，所述分層貝葉斯模型的目標(biāo)函數(shù)中包括表示微生物之間的關(guān)聯(lián)的精確度矩 陣和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣；
[0021 ] 相應(yīng)地，所述對所述目標(biāo)函數(shù)進行優(yōu)化，包括：
[0022]利用graphical lasso算法對表示微生物之間的關(guān)聯(lián)的精確度矩陣進行迭代優(yōu) 化；
[0023]利用proximal算法對表示微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣進行迭代 優(yōu)化。
[0024]其中，在所述獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組測 序樣本對應(yīng)的環(huán)境因素測量值之前，所述方法還包括：
[0025]對所述宏基因組測序樣本進行預(yù)處理。
[0026]根據(jù)本發(fā)明的第二方面，提供了一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測裝置，該裝置包括：
[0027] 獲取單元，用于獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組 測序樣本對應(yīng)的環(huán)境因素測量值；
[0028] 模型建立單元，用于根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組 測序樣本中各種微生物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型，所述分 層貝葉斯模型用于對數(shù)據(jù)產(chǎn)生過程、微生物之間的關(guān)聯(lián)以及微生物與環(huán)境因素測量值之間 的關(guān)聯(lián)進行描述；
[0029]模型學(xué)習(xí)單元，用于采用最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確 定所述分層貝葉斯模型的目標(biāo)函數(shù)；
[0030] 優(yōu)化單元，用于對所述目標(biāo)函數(shù)進行優(yōu)化；
[0031] 預(yù)測單元，用于采用目標(biāo)函數(shù)優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和 微生物與環(huán)境因素測量值之間的關(guān)聯(lián)。
[0032] 其中，所述模型建立單元，具體用于采用Dirichlet-Multinomial共輒分布模擬測 序過程，對數(shù)據(jù)產(chǎn)生過程中的組成成分偏差和數(shù)據(jù)方差進行第一建模;采用Lognormal分布 模擬微生物的豐度變化過程，對微生物之間的關(guān)聯(lián)以及微生物與環(huán)境因素測量值之間的關(guān) 聯(lián)進行第二建模;根據(jù)所述第一建模和第二建模的結(jié)果，建立基于Lognormal-Dirichlet-Multinomial的分層貝葉斯模型。
[0033] 其中，所述模型學(xué)習(xí)單元，具體用于采用增加一范數(shù)懲罰項的最大后驗估計算法 對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分層貝葉斯模型的目標(biāo)函數(shù)。
[0034] 其中，所述分層貝葉斯模型的目標(biāo)函數(shù)中包括表示微生物之間的關(guān)聯(lián)的精確度矩 陣和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣；
[0035] 相應(yīng)地，所述優(yōu)化單元，包括：
[0036 ]第一優(yōu)化模塊，用于利用graph i c a I I a s s 〇算法對表示微生物之間的關(guān)聯(lián)的精確 度矩陣進行迭代優(yōu)化；
[0037]第二優(yōu)化模塊，用于利用proximal算法對表示微生物與環(huán)境因素測量值之間的關(guān) 聯(lián)的矩陣進行迭代優(yōu)化。
[0038] 其中，所述裝置還包括：
[0039]預(yù)處理單元，用于在所述獲取宏基因組測序樣本中各種微生物的豐度，以及每一 宏基因組測序樣本對應(yīng)的環(huán)境因素測量值之前，對所述宏基因組測序樣本進行預(yù)處理。
[0040] 本發(fā)明提供的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法及裝置，根據(jù)宏基因組測序數(shù)據(jù)的產(chǎn)生過 程以及每一宏基因組測序樣本中各種微生物的豐度和對應(yīng)的環(huán)境因素測量值，建立分層貝 葉斯模型，并對該模型進行學(xué)習(xí)和優(yōu)化，通過采用優(yōu)化后的分層貝葉斯模型預(yù)測微生物之 間的關(guān)聯(lián)和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)，以解決由于組成成分偏差和測序數(shù)據(jù)自 身的方差所導(dǎo)致的關(guān)聯(lián)推斷不準(zhǔn)確的問題，同時考慮微生物和環(huán)境因素的影響，顯著提高 在微生物關(guān)聯(lián)和微生物與環(huán)境因素關(guān)聯(lián)的預(yù)測任務(wù)中的準(zhǔn)確性和實用性。
【附圖說明】
[0041] 通過閱讀下文優(yōu)選實施方式的詳細(xì)描述，各種其他的優(yōu)點和益處對于本領(lǐng)域普通 技術(shù)人員將變得清楚明了。附圖僅用于示出優(yōu)選實施方式的目的，而并不認(rèn)為是對本發(fā)明 的限制。而且在整個附圖中，用相同的參考符號表示相同的部件。在附圖中：
[0042] 圖1為本發(fā)明一個實施例提出的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法的流程圖；
[0043] 圖2為本發(fā)明另一實施例提出的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法的流程圖；
[0044] 圖3為本發(fā)明一個實施例提出的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測裝置的結(jié)構(gòu)示意圖；
[0045] 圖4為本發(fā)明另一實施例提出的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測裝置的結(jié)構(gòu)示意圖。
【具體實施方式】
[0046] 下面詳細(xì)描述本發(fā)明的實施例，所述實施例的示例在附圖中示出，其中自始至終 相同或類似的標(biāo)號表示相同或類似的元件或具有相同或類似功能的元件。下面通過參考附 圖描述的實施例是示例性的，僅用于解釋本發(fā)明，而不能解釋為對本發(fā)明的限制。
[0047] 圖1示出了本發(fā)明實施例的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法的流程圖。參照圖1，本發(fā)明 實施例提出的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法，具體包括以下步驟：
[0048] 步驟Sll、獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組測序樣 本對應(yīng)的環(huán)境因素測量值。
[0049]在實際應(yīng)用中，通過對宏基因組數(shù)據(jù)進行處理，獲取樣本中各種微生物(OTU)的豐 度和環(huán)境因素的值，同時去掉異常的樣本和特征，便于模型的下一步的計算。對于宏基因組 16S測序數(shù)據(jù)，需要對樣本進行質(zhì)量控制，reads聚類和注釋，然后才能得到OTU的豐度;環(huán)境 因素需要在取樣的時候就已經(jīng)測定。
[0050] 步驟Sl 2、根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組測序樣本中 各種微生物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型，所述分層貝葉斯模 型用于對數(shù)據(jù)產(chǎn)生過程、微生物之間的關(guān)聯(lián)以及微生物與環(huán)境因素測量值之間的關(guān)聯(lián)進行 描述。
[0051] 步驟S13、采用最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分層 貝葉斯模型的目標(biāo)函數(shù)；
[0052]步驟S14、對所述目標(biāo)函數(shù)進行優(yōu)化；
[0053]步驟S15、采用目標(biāo)函數(shù)優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和微生 物與環(huán)境因素測量值之間的關(guān)聯(lián)。
[0054]本發(fā)明實施例提供的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法，根據(jù)宏基因組測序數(shù)據(jù)的產(chǎn)生過 程以及每一宏基因組測序樣本中各種微生物的豐度和對應(yīng)的環(huán)境因素測量值，建立分層貝 葉斯模型，并對該模型進行學(xué)習(xí)和優(yōu)化，通過采用優(yōu)化后的分層貝葉斯模型預(yù)測微生物之 間的關(guān)聯(lián)和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)，以解決由于組成成分偏差和測序數(shù)據(jù)自 身的方差所導(dǎo)致的關(guān)聯(lián)推斷不準(zhǔn)確的問題，同時考慮微生物和環(huán)境因素的影響，顯著提高 在微生物關(guān)聯(lián)和微生物與環(huán)境因素關(guān)聯(lián)的預(yù)測任務(wù)中的準(zhǔn)確性和實用性。
[0055]在本發(fā)明的一個可選實施例中，步驟S12,具體包括以下附圖中未示出的步驟： [0056] 步驟S121、采用Dirichlet-Multinomial共輒分布模擬測序過程，對數(shù)據(jù)產(chǎn)生過程 中的組成成分偏差和數(shù)據(jù)方差進行第一建模；
[0057]步驟S122、采用Lognormal分布模擬微生物的豐度變化過程，對微生物之間的關(guān)聯(lián) 以及微生物與環(huán)境因素測量值之間的關(guān)聯(lián)進行第二建模；
[0058] 步驟S123、根據(jù)所述第一建模和第二建模的結(jié)果，建立基于Lognormal-Di;richlet-Multinomial的分層貝葉斯模型。
[0059] 在實際應(yīng)用中，根據(jù)宏基因組測序數(shù)據(jù)的特點，建立分層貝葉斯模型，對數(shù)據(jù)產(chǎn)生 過程、微生物之間的關(guān)聯(lián)和微生物與環(huán)境因素之間的關(guān)聯(lián)進行描述;為了考慮測序數(shù)據(jù)樣 本之間reads總數(shù)不同，簡單的歸一化會引入組成成分偏差；以及測序數(shù)據(jù)由于一系列處理 步驟所產(chǎn)生的誤差的問題，本發(fā)明用Dirichlet-Multinomial共輒分布進行模擬。假設(shè)測序 得到的微生物的reads數(shù)服從多項式分布，這些reads數(shù)受到樣本的大小和該OTU所占相對 比例的影響，進一步假設(shè)其相對比例服從狄利克雷分布。OTU的相對豐度進一步由其環(huán)境中 對應(yīng)的微生物的絕度豐度決定，而微生物的絕對豐度假設(shè)受到兩方面因素的影響，一方面 是微生物之間的關(guān)聯(lián)，另一方面是微生物與環(huán)境因素之間的關(guān)聯(lián)，本發(fā)明對這一過程用 Lognormal進行建模。最后，微生物之間的關(guān)聯(lián)和微生物與環(huán)境因素的分別參數(shù)化為兩個矩 陣，這兩個矩陣對應(yīng)Lognormal分布的參數(shù)。
[0060] 下面通過一個具體實施例，對本發(fā)明技術(shù)方案中根據(jù)每一宏基因組測序樣本的數(shù) 據(jù)產(chǎn)生過程、每一宏基因組測序樣本中各種微生物的豐度以及對應(yīng)的環(huán)境因素測量值，建 立分層貝葉斯模型的實現(xiàn)方式進行清楚的說明。
[0061 ] 假設(shè)Xi是一個P維的向量，由一個樣本中P個OTU測序得到的reads數(shù)組成，mi是一個 Q維的向量，由Q個環(huán)境因素的值構(gòu)成;hi和Xi相對應(yīng)，表示的是P個OTU對應(yīng)的在樣本中的相 對豐度;Q i表示的是P個OTUs對應(yīng)的在真實環(huán)境中的絕對豐度，Z表示隱變量，絕對豐度又收 到兩方面的影響:一方面是微生物之間的關(guān)聯(lián)，用Θ精確度矩陣表示(P X P矩陣），另一方面 是環(huán)境因素的影響，用矩陣B表示(QXP矩陣）。此是對應(yīng)P個OTU豐度的基向量。
[0062] Zi~Gaussian(Bt), Θ-1)
[0063] Cii = exp (BTmi+Zi)
[0064] hi ~Dirichlet(ai)
[0065] Xi ~Multinomial(hi)
[0066] 因為微生物的reads ^是通過測序得到，測序的過程中的PCR可以用多項式分布 來建模，即在樣本總的reads數(shù)確定時，每種OTU的數(shù)量和其相對比例有關(guān)，相當(dāng)于按其相對 比例進行采樣：
[0067]
[0068] 其中
§第i個樣本中的總reads數(shù)。hi對應(yīng)P個OTUs的相對分
微生物的相對豐度實際上是根據(jù)其絕對豐度，即在群落中的絕對數(shù)量計算得到的，本發(fā)明 用Dirchlet分布來建模相對分度hi和絕對豐度Cti之間的關(guān)系，
[0069] )
[0070] ]馬函數(shù)。利用Dirichlet-Multinomial共輒分布的性
質(zhì)，
[0071
[0072] 可得OTU-j在第i個樣本中的方差為Var(Xij) = S(Xi) · C · rij · (1-rij)，和第j個 OTU的協(xié)方差為Cov(Xij，Xik)=_s(Xi) · C · rij · rik，其中c ，.:= 和rik = alk/s(ai)。可以看到利用該共輒分布，微生物的數(shù)量之間有了一個負(fù)的相關(guān)性，這個正好 與組成成分偏差對應(yīng);同時，其方差和協(xié)方差均與絕對豐度ai有關(guān)，rij正好是由絕對豐度決 定的相對豐度的期望值，與測序數(shù)據(jù)較為吻合。
[0073]進一步的，本發(fā)明假設(shè)微生物的絕對豐度服從Lognormal分布，這個分布被生物學(xué) 家廣泛采用。同時，為了考慮微生物和微生物之間的關(guān)聯(lián)和環(huán)境因素對微生物的影響，我們 I ^4V /?￥Τ K 1 固 .
[0074]
[0075] 其中，均值yi = BV+Bo。這樣，精確度矩陣Θ就用來表征微生物之間的關(guān)聯(lián)，矩陣B 用來表示微生物和環(huán)境因素之間的關(guān)聯(lián)。
[0076] 在本發(fā)明的一個可選實施例中，步驟S13,具體包括以下附圖中未示出的步驟:采 用增加一范數(shù)懲罰項的最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分層 貝葉斯模型的目標(biāo)函數(shù)。
[0077] 在實際應(yīng)用中，在Lognormal-Dirichlet-Multinomial分層貝葉斯模型的基礎(chǔ)上， 對微生物之間的關(guān)聯(lián)和微生物與環(huán)境因素之間的關(guān)聯(lián)引入一范數(shù)的稀疏約束，利用最大后 驗進行估計。對要推斷的關(guān)聯(lián)進行稀疏約束，一方面是因為微生物的數(shù)量遠多于一個群落 可以獲得的樣本數(shù)，加入稀疏約束可以處理過擬合的問題；另一方面是因為測序數(shù)據(jù)噪聲 較大，通過稀疏項可以保留最顯著的關(guān)聯(lián)，提高準(zhǔn)確性。由于本發(fā)明的分層貝葉斯模型中含 有隱變量，在利用最大后驗進行參數(shù)估計的過程中，隱變量也一同進行估計，這樣降低了優(yōu) 化的復(fù)雜性。
[0078] 下面通過一個具體實施例，對本發(fā)明技術(shù)方案中采用最大后驗估計算法對所述分 層貝葉斯模型進行學(xué)習(xí)，確定所述分層貝葉斯模型的目標(biāo)函數(shù)的實現(xiàn)方式進行清楚的說 明。
[0079] 采用增加一范數(shù)懲罰項的最大后驗估計方法學(xué)習(xí)該分層貝葉斯模型。由于本發(fā)明 中的分層模型含有隱變量，直接估計較為困難，所以采用最大后驗，帶著隱變量一起估計。 假設(shè)宏基因組測序OTU的數(shù)據(jù)為矩陣Χ(Ν*Ρ矩陣，N個樣本，P個0TU)，環(huán)境數(shù)據(jù)為矩陣M(N*Q 矩陣，N個樣本，Q個環(huán)境因素），則對于隱變量Z的最大后驗估計為：
[0080] P(Z|X,M,B,Bo,Q)^P(X,Z|B,Bo,Q,M)^P(X|a)P(a|Z,B,Bo,M)P(Z|Bo,Q)
[0081] 其中P(X|a)可以用等式（3)計算得到，并且I3)，P(Zl | Bo，Q)是高 斯分布。

[0082] 采用一范數(shù)一方面是因為模型參數(shù)較多，樣本較少，需要利用一范數(shù)懲罰項防止 過擬合;另一方面是便于做特征選擇，保留出最顯著的關(guān)聯(lián)供以后的深入分析。最終，本發(fā) 明的曰標(biāo)函撒為-
[008：
[0084]其中，
[0086] 在本發(fā)明的一個可選實施例中，所述分層貝葉斯模型的目標(biāo)函數(shù)中包括表示微生 物之間的關(guān)聯(lián)的精確度矩陣和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣。
[0087] 在本發(fā)明實施例中，步驟S14,具體包括以下附圖中未示出的步驟：
[0088]步驟S141、利用graphical lasso算法對表示微生物之間的關(guān)聯(lián)的精確度矩陣進 行迭代優(yōu)化；
[0089]步驟S142、利用proximal算法對表示微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣 進行迭代優(yōu)化。
[0090] 具體的，本發(fā)明實施例對目標(biāo)函數(shù)利用proximal方法和graphical lasso (Friedman et.al.，2008)方法進行迭代優(yōu)化。對于表示微生物之間關(guān)聯(lián)的精確度矩陣，利 用graphical lasso進行優(yōu)化，這種方法快速有效;對于表示微生物與環(huán)境因素之間關(guān)聯(lián)的 矩陣，我們用proximal方法進行優(yōu)化，每次迭代都對目標(biāo)函數(shù)做二階近似，然后通過 coordinate descent坐標(biāo)下降法進行優(yōu)化。
[0091]下面通過一個具體實施例，對本發(fā)明技術(shù)方案中對所述目標(biāo)函數(shù)進行優(yōu)化的實現(xiàn) 方式進行清楚的說明。
[0092] 本發(fā)明對目標(biāo)函數(shù)(7)利用proximal方法和graphical lasso方法進行迭代優(yōu)化。 由于模型中未知參數(shù)包括:Z，B，Q和Bo,本發(fā)明采用分塊迭代的方法，每次優(yōu)化一個參數(shù)，然 后通過不斷交替迭代最終收斂。
[0093] 對于隱變量Z，對目標(biāo)函數(shù)導(dǎo)數(shù)為：
[0094]
[0095]可以利用L-BFGS擬牛頓法進行優(yōu)化。
[0096] 對于矩陣B，利用proximal擬牛頓方法進行迭代。對于矩陣B，導(dǎo)數(shù)如下：
[0097；
[0098]然后利用一階導(dǎo)數(shù)近似Hessian矩陣，從而得到二階近似目標(biāo)函數(shù)，加上一范數(shù)的 約束。對于向量Bo，更新規(guī)則為：
[0099]
[0100] 對于矩陣Q，其目標(biāo)函數(shù)為：
[0101 ] nun- log | Q j +ir(SQ] + Λ, ]j Q |], (11)
[0102] 該目標(biāo)函數(shù)可以通過graphical lasso算法有效迭代求解。
[0103] 綜上，基于Lognormal-Dir ichle t-Multinomial分層貝葉斯模型的核心學(xué)習(xí)算法 如下：
[0106] 本發(fā)明實施例估計出的兩個參數(shù):矩陣Q和B，分別用于解釋和預(yù)測微生物之間的 關(guān)聯(lián)和微生物與環(huán)境因素之間的關(guān)聯(lián)。矩陣中的元素叫表示了 〇TU-i和OTU-j之間的關(guān)聯(lián)， 如果Qu = 〇,則〇TU-i和j是條件獨立的，否則是條件相關(guān)，該關(guān)聯(lián)權(quán)重>
對于 矩陣B，B^表示OTU-j和環(huán)境因素 i之間的關(guān)聯(lián)，其權(quán)值為在估計出關(guān)聯(lián)后，可以利用關(guān) 聯(lián)的正負(fù)性和絕對值的大小，結(jié)合OTU注釋的微生物物種信息，推薦二者之間真實的交互關(guān) 系，如競爭，共生和捕食等。
[0107] 圖2示出了本發(fā)明實施例的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法的流程圖。參照圖2,本發(fā)明 實施例提出的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法，在所述獲取宏基因組測序樣本中各種微生物的豐 度，以及每一宏基因組測序樣本對應(yīng)的環(huán)境因素測量值之前，所述方法還包括步驟S10:
[0108] 步驟SI 0、對所述宏基因組測序樣本進行預(yù)處理。
[0109] 在具體應(yīng)用中，由于實際采樣和測序過程的誤差，會導(dǎo)致某些環(huán)境因素的值未測 量或明顯與實際不符，某些OTU的reads數(shù)量也會出現(xiàn)較大的波動，樣本之間總reads數(shù)也會 出現(xiàn)很大的差別，這些因素都需要考慮進來，對得到的數(shù)據(jù)進行過濾等預(yù)處理。這樣得到的 數(shù)據(jù)作為本發(fā)明模型的輸入。
[0110]對于方法實施例，為了簡單描述，故將其都表述為一系列的動作組合，但是本領(lǐng)域 技術(shù)人員應(yīng)該知悉，本發(fā)明實施例并不受所描述的動作順序的限制，因為依據(jù)本發(fā)明實施 例，某些步驟可以采用其他順序或者同時進行。其次，本領(lǐng)域技術(shù)人員也應(yīng)該知悉，說明書 中所描述的實施例均屬于優(yōu)選實施例，所涉及的動作并不一定是本發(fā)明實施例所必須的。
[0111]圖3示出了本發(fā)明實施例的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測裝置的結(jié)構(gòu)示意圖。參照圖3,本 發(fā)明實施例提供的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測裝置，包括獲取單元301、模型建立單元302、模型學(xué) 習(xí)單元303、優(yōu)化單元304以及預(yù)測單元305,其中：獲取單元301用于獲取宏基因組測序樣本 中各種微生物的豐度，以及每一宏基因組測序樣本對應(yīng)的環(huán)境因素測量值;模型建立單元 302用于根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組測序樣本中各種微生 物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型，所述分層貝葉斯模型用于對 數(shù)據(jù)產(chǎn)生過程、微生物之間的關(guān)聯(lián)以及微生物與環(huán)境因素測量值之間的關(guān)聯(lián)進行描述;模 型學(xué)習(xí)單元303用于采用最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分 層貝葉斯模型的目標(biāo)函數(shù);優(yōu)化單元304用于對所述目標(biāo)函數(shù)進行優(yōu)化;預(yù)測單元305用于 采用目標(biāo)函數(shù)優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和微生物與環(huán)境因素測量 值之間的關(guān)聯(lián)。
[0112]本發(fā)明提供的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測裝置，模型建立單元302根據(jù)宏基因組測序數(shù) 據(jù)的產(chǎn)生過程以及每一宏基因組測序樣本中各種微生物的豐度和對應(yīng)的環(huán)境因素測量值， 建立分層貝葉斯模型，并通過模型學(xué)習(xí)單元303和優(yōu)化單元304對該模型進行學(xué)習(xí)和優(yōu)化， 以供預(yù)測單元305采用優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和微生物與環(huán)境因 素測量值之間的關(guān)聯(lián)，以解決由于組成成分偏差和測序數(shù)據(jù)自身的方差所導(dǎo)致的關(guān)聯(lián)推斷 不準(zhǔn)確的問題，同時考慮微生物和環(huán)境因素的影響，顯著提高在微生物關(guān)聯(lián)和微生物與環(huán) 境因素關(guān)聯(lián)的預(yù)測任務(wù)中的準(zhǔn)確性和實用性。
[0113]在本發(fā)明的一個可選實施例中，所述模型建立單元302,具體用于采用Dirichlet-Multinomial共輒分布模擬測序過程，對數(shù)據(jù)產(chǎn)生過程中的組成成分偏差和數(shù)據(jù)方差進行 第一建模;采用Lognormal分布模擬微生物的豐度變化過程，對微生物之間的關(guān)聯(lián)以及微生 物與環(huán)境因素測量值之間的關(guān)聯(lián)進行第二建模;根據(jù)所述第一建模和第二建模的結(jié)果，建 立基于Lognormal-Dir i chi et-Mult inomial的分層貝葉斯模型。
[0114] 在本發(fā)明的一個可選實施例中，所述模型學(xué)習(xí)單元303,具體用于采用增加一范數(shù) 懲罰項的最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分層貝葉斯模型的 目標(biāo)函數(shù)。
[0115] 在本發(fā)明的一個可選實施例中，所述分層貝葉斯模型的目標(biāo)函數(shù)中包括表示微生 物之間的關(guān)聯(lián)的精確度矩陣和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣；
[0116]相應(yīng)地，所述優(yōu)化單元304，具體包括第一優(yōu)化模塊和第二優(yōu)化模塊，其中，第一優(yōu) 化模塊，用于利用graphical lasso算法對表示微生物之間的關(guān)聯(lián)的精確度矩陣進行迭代 優(yōu)化;第二優(yōu)化模塊，用于利用proximal算法對表示微生物與環(huán)境因素測量值之間的關(guān)聯(lián) 的矩陣進行迭代優(yōu)化。
[0117] 在本發(fā)明的一個可選實施例中，如圖4所示，所述裝置還包括預(yù)處理單元300，所述 預(yù)處理單元300,用于在所述獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因 組測序樣本對應(yīng)的環(huán)境因素測量值之前，對所述宏基因組測序樣本進行預(yù)處理。
[0118] 對于裝置實施例而言，由于其與方法實施例基本相似，所以描述的比較簡單，相關(guān) 之處參見方法實施例的部分說明即可。
[0119] 綜上所述，本發(fā)明實施例提供的微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法及裝置，根據(jù)宏基因組 測序數(shù)據(jù)的產(chǎn)生過程以及每一宏基因組測序樣本中各種微生物的豐度和對應(yīng)的環(huán)境因素 測量值，建立分層貝葉斯模型，并對該模型進行學(xué)習(xí)和優(yōu)化，通過采用優(yōu)化后的分層貝葉斯 模型預(yù)測微生物之間的關(guān)聯(lián)和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)，以解決由于組成成分 偏差和測序數(shù)據(jù)自身的方差所導(dǎo)致的關(guān)聯(lián)推斷不準(zhǔn)確的問題，同時考慮微生物和環(huán)境因素 的影響，顯著提高在微生物關(guān)聯(lián)和微生物與環(huán)境因素關(guān)聯(lián)的預(yù)測任務(wù)中的準(zhǔn)確性和實用 性。
[0120] 通過以上的實施方式的描述，本領(lǐng)域的技術(shù)人員可以清楚地了解到本發(fā)明可以通 過硬件實現(xiàn)，也可以借助軟件加必要的通用硬件平臺的方式來實現(xiàn)?；谶@樣的理解，本發(fā) 明的技術(shù)方案可以以軟件產(chǎn)品的形式體現(xiàn)出來，該軟件產(chǎn)品可以存儲在一個非易失性存儲 介質(zhì)(可以是CD-R0M，U盤，移動硬盤等）中，包括若干指令用以使得一臺計算機設(shè)備(可以是 個人計算機，服務(wù)器，或者網(wǎng)絡(luò)設(shè)備等)執(zhí)行本發(fā)明各個實施例所述的方法。
[0121] 本領(lǐng)域技術(shù)人員可以理解附圖只是一個優(yōu)選實施例的示意圖，附圖中的模塊或流 程并不一定是實施本發(fā)明所必須的。
[0122] 本領(lǐng)域技術(shù)人員可以理解實施例中的裝置中的模塊可以按照實施例描述進行分 布于實施例的裝置中，也可以進行相應(yīng)變化位于不同于本實施例的一個或多個裝置中。上 述實施例的單元可以合并為一個單元，也可以進一步拆分成多個子模塊。
[0123] 以上所述僅是本發(fā)明的部分實施方式，應(yīng)當(dāng)指出，對于本技術(shù)領(lǐng)域的普通技術(shù)人 員來說，在不脫離本發(fā)明原理的前提下，還可以做出若干改進和潤飾，這些改進和潤飾也應(yīng) 視為本發(fā)明的保護范圍。
【主權(quán)項】
1. 一種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測方法，其特征在于，該方法包括： 獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組測序樣本對應(yīng)的環(huán)境 因素測量值； 根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組測序樣本中各種微生物的 豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型，所述分層貝葉斯模型用于對數(shù)據(jù) 產(chǎn)生過程、微生物之間的關(guān)聯(lián)以及微生物與環(huán)境因素測量值之間的關(guān)聯(lián)進行描述； 采用最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分層貝葉斯模型的 目標(biāo)函數(shù)； 對所述目標(biāo)函數(shù)進行優(yōu)化； 采用目標(biāo)函數(shù)優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和微生物與環(huán)境因素 測量值之間的關(guān)聯(lián)。2. 根據(jù)權(quán)利要求1所述的方法，其特征在于，所述根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn) 生過程、每一宏基因組測序樣本中各種微生物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分 層貝葉斯模型，包括： 采用0;^;[(311161：-]/[1111:;[11〇1]11&1共輒分布模擬測序過程，對數(shù)據(jù)產(chǎn)生過程中的組成成分 偏差和數(shù)據(jù)方差進行第一建模； 采用Lognormal分布模擬微生物的豐度變化過程，對微生物之間的關(guān)聯(lián)以及微生物與 環(huán)境因素測量值之間的關(guān)聯(lián)進行第二建模； 根據(jù)所述第一建模和第二建模的結(jié)果，建立基于Lognormal-Dirichlet-Multinomial 的分層貝葉斯模型。3. 根據(jù)權(quán)利要求1所述的方法，其特征在于，所述采用最大后驗估計算法對所述分層貝 葉斯模型進行學(xué)習(xí)，確定所述分層貝葉斯模型的目標(biāo)函數(shù)，包括： 采用增加一范數(shù)懲罰項的最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所 述分層貝葉斯模型的目標(biāo)函數(shù)。4. 根據(jù)權(quán)利要求1所述的方法，其特征在于，所述分層貝葉斯模型的目標(biāo)函數(shù)中包括表 示微生物之間的關(guān)聯(lián)的精確度矩陣和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣； 相應(yīng)地，所述對所述目標(biāo)函數(shù)進行優(yōu)化，包括： 利用graphical lasso算法對表示微生物之間的關(guān)聯(lián)的精確度矩陣進行迭代優(yōu)化； 利用proximal算法對表示微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣進行迭代優(yōu)化。5. 根據(jù)權(quán)利要求1-4任一項所述的方法，其特征在于，在所述獲取宏基因組測序樣本中 各種微生物的豐度，以及每一宏基因組測序樣本對應(yīng)的環(huán)境因素測量值之前，所述方法還 包括： 對所述宏基因組測序樣本進行預(yù)處理。6. -種微生物關(guān)聯(lián)網(wǎng)絡(luò)預(yù)測裝置，其特征在于，該裝置包括： 獲取單元，用于獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組測序 樣本對應(yīng)的環(huán)境因素測量值； 模型建立單元，用于根據(jù)每一宏基因組測序樣本的數(shù)據(jù)產(chǎn)生過程、每一宏基因組測序 樣本中各種微生物的豐度以及對應(yīng)的環(huán)境因素測量值，建立分層貝葉斯模型，所述分層貝 葉斯模型用于對數(shù)據(jù)產(chǎn)生過程、微生物之間的關(guān)聯(lián)以及微生物與環(huán)境因素測量值之間的關(guān) 聯(lián)進行描述； 模型學(xué)習(xí)單元，用于采用最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所 述分層貝葉斯模型的目標(biāo)函數(shù)； 優(yōu)化單元，用于對所述目標(biāo)函數(shù)進行優(yōu)化； 預(yù)測單元，用于采用目標(biāo)函數(shù)優(yōu)化后的分層貝葉斯模型預(yù)測微生物之間的關(guān)聯(lián)和微生 物與環(huán)境因素測量值之間的關(guān)聯(lián)。7. 根據(jù)權(quán)利要求6所述的裝置，其特征在于，所述模型建立單元，具體用于采用 Dirichlet-Multinomial共輒分布模擬測序過程，對數(shù)據(jù)產(chǎn)生過程中的組成成分偏差和數(shù) 據(jù)方差進行第一建模;采用Lognormal分布模擬微生物的豐度變化過程，對微生物之間的關(guān) 聯(lián)以及微生物與環(huán)境因素測量值之間的關(guān)聯(lián)進行第二建模;根據(jù)所述第一建模和第二建模 的結(jié)果，建立基于Lognormal-Dirichlet-Multinomial的分層貝葉斯模型。8. 根據(jù)權(quán)利要求6所述的裝置，其特征在于，所述模型學(xué)習(xí)單元，具體用于采用增加一 范數(shù)懲罰項的最大后驗估計算法對所述分層貝葉斯模型進行學(xué)習(xí)，確定所述分層貝葉斯模 型的目標(biāo)函數(shù)。9. 根據(jù)權(quán)利要求6所述的裝置，其特征在于，所述分層貝葉斯模型的目標(biāo)函數(shù)中包括表 示微生物之間的關(guān)聯(lián)的精確度矩陣和微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的矩陣； 相應(yīng)地，所述優(yōu)化單元，包括： 第一優(yōu)化模塊，用于利用graphical lasso算法對表示微生物之間的關(guān)聯(lián)的精確度矩 陣進行迭代優(yōu)化； 第二優(yōu)化模塊，用于利用proximal算法對表示微生物與環(huán)境因素測量值之間的關(guān)聯(lián)的 矩陣進行迭代優(yōu)化。10. 根據(jù)權(quán)利要求6-9任一項所述的裝置，其特征在于，所述裝置還包括:預(yù)處理單元， 用于在所述獲取宏基因組測序樣本中各種微生物的豐度，以及每一宏基因組測序樣本對應(yīng) 的環(huán)境因素測量值之前，對所述宏基因組測序樣本進行預(yù)處理。
【文檔編號】G06F19/24GK105938524SQ201610266864
【公開日】2016年9月14日
【申請日】2016年4月26日
【發(fā)明人】陳挺, 陳寧, 楊煜清
【申請人】清華大學(xué)