多組學(xué)數(shù)據(jù)分析中關(guān)聯(lián)網(wǎng)絡(luò)融合ANF算法及應(yīng)用研究*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計學(xué)教研室(150081) 徐臻旖 王 策 侯 艷 李 康

【提 要】 目的 引入關(guān)聯(lián)網(wǎng)絡(luò)融合(affinity network fusion，ANF)方法對多組學(xué)數(shù)據(jù)進(jìn)行整合分析，并應(yīng)用于腫瘤分子分型研究。方法 模擬產(chǎn)生兩組學(xué)數(shù)據(jù)，改變總體差異大小等情況比較多種多組學(xué)整合方法的效果。實例分析中目標(biāo)人群選擇TCGA數(shù)據(jù)庫中對鉑類藥物敏感并擁有mRNA和甲基化兩個組學(xué)的卵巢癌患者，目標(biāo)基因是TCGA數(shù)據(jù)庫和ImmPort數(shù)據(jù)庫中共有基因，目標(biāo)甲基化位點是目標(biāo)基因?qū)?yīng)的所有甲基化位點。使用ANF、SNF、K-means、系統(tǒng)聚類和iCluster五種方法比較聚類效果。結(jié)果 模擬實驗提示存在總體差異的兩亞型間差異較小時ANF方法的效果明顯優(yōu)于其他方法。實例分析結(jié)果表明，通過ANF方法進(jìn)行多組學(xué)數(shù)據(jù)整合得到的分子分型較單組學(xué)得到的分子分型有更好的生物學(xué)意義且多組學(xué)聚類效果優(yōu)于其他方法。結(jié)論 ANF方法可以應(yīng)用于多組學(xué)數(shù)據(jù)整合分析，具有實際應(yīng)用意義。

多組學(xué)指基因組、轉(zhuǎn)錄組、蛋白質(zhì)組及代謝組等，利用多組學(xué)進(jìn)行整合分析能提供更多的信息。傳統(tǒng)的分析方法是對各組學(xué)分別進(jìn)行分析，然后再根據(jù)生物學(xué)知識對其進(jìn)行解釋。這種分析方法的不足之處在于不能充分利用數(shù)據(jù)提供的信息，從根本上建立各組學(xué)之間的關(guān)系。根據(jù)不同目的可將多組學(xué)整合方法分為不同類型，如基于奇異值分解發(fā)現(xiàn)組學(xué)間調(diào)控關(guān)系的JIVE(joint and individual variation explained)方法[1]，基于聯(lián)合潛在變量模型的iCluster方法以及基于樣本間網(wǎng)絡(luò)結(jié)構(gòu)得到分子亞型的相似網(wǎng)絡(luò)融合(similarity network fusion，SNF)方法等[2-3]。本文引進(jìn)關(guān)聯(lián)網(wǎng)絡(luò)融合法(affinity network fusion，ANF)，與SNF方法相比計算過程更簡潔，可利用先驗信息對不同組學(xué)進(jìn)行加權(quán)。ANF將不同組學(xué)的數(shù)據(jù)矩陣進(jìn)行轉(zhuǎn)化再整合分析，為多組學(xué)研究提供有價值的信息[4]。

原理和方法

ANF是一種利用多組學(xué)數(shù)據(jù)構(gòu)建矩陣的方法，即通過隨機(jī)游走的方式將分別構(gòu)建的組學(xué)樣本間關(guān)聯(lián)矩陣進(jìn)一步整合[4]，對最終得到的包含多組學(xué)信息的關(guān)聯(lián)矩陣進(jìn)行譜聚類分析。

具體步驟如下：

1.構(gòu)建關(guān)聯(lián)矩陣

(1)

(2)

(3)

2.不同組學(xué)關(guān)聯(lián)矩陣融合

(4)

(5)

對得到的關(guān)聯(lián)矩陣W進(jìn)行多分類譜聚類分析[5]。多分類譜聚類由Stella等人提出，與標(biāo)準(zhǔn)譜聚類不同在于該方法使用奇異值分解進(jìn)行計算，聚類數(shù)目通過啟發(fā)式特征間隙法(eigengap heuristic)進(jìn)行確定，假設(shè)求得m個特征值λ1，…，λm，存在i∈[1，m]，使得λi與λi+1的差值較大，則認(rèn)為聚類類別數(shù)為i類。

3.聚類效果評價

聚類效果的衡量指標(biāo)通常分為兩大類：第一類為有監(jiān)督方法，已知真實標(biāo)簽，用一定的度量評判聚類結(jié)果與真實標(biāo)簽的符合程度，例如標(biāo)準(zhǔn)化互信息(NMI)[0，1]，NMI的值越接近1表明得到聚類結(jié)果與真實情況越吻合；第二類為無監(jiān)督方法，根據(jù)是否滿足類內(nèi)聚集程度高與類間離散程度大對聚類效果進(jìn)行評價，例如輪廓系數(shù)(silhouette coefficient)，其值越接近1表明聚類越合理，越接近-1表示可能存在錯分情況；還可通過比較生存分析的預(yù)后情況評價聚類效果。

模擬實驗

可以看到，隨著差異變量比例和不同亞組差異變量均值差異增加，NMI逐漸增加，ANF在不同亞型總體差異較小時效果優(yōu)于其他四種方法，如均數(shù)差值為4，變量數(shù)為200，差異變量為15%時，ANF的NMI值為0.750，明顯大于其他聚類方法的NMI值。當(dāng)不同亞型總體差異逐漸增加達(dá)到一定程度時，除系統(tǒng)聚類外以上方法的聚類效果均較好。

實例分析

本研究對使用鉑類藥物敏感的卵巢癌患者進(jìn)行免疫亞型分析。選擇TCGA數(shù)據(jù)庫中同時擁有mRNA、甲基化數(shù)據(jù)以及患者臨床信息的卵巢癌樣本，使用ImmPort數(shù)據(jù)庫確定與免疫相關(guān)的mRNA?；熋舾谢颊叨x為使用大于等于六個周期的鉑類藥物且化療后無進(jìn)展生存期為六個月以上的患者；將化療耐藥患者定義為第一次使用鉑類藥物化療耐藥以及在治療周期完成后六個月內(nèi)復(fù)發(fā)的患者[7]。

對數(shù)據(jù)進(jìn)行預(yù)處理，將樣本ID匹配，要求樣本同時具有mRNA數(shù)據(jù)和甲基化數(shù)據(jù)，選擇TCGA與ImmPort數(shù)據(jù)庫中共有的mRNA，同時確定mRNA所對應(yīng)的甲基化位點。刪除重復(fù)變量，將mRNA中超過70%為0的變量剔除，將甲基化位點中超過70%缺失的變量剔除，并用最近鄰KNN算法填補(bǔ)剩余缺失值。最終鉑類藥物化療敏感人群136人，1211個mRNA以及2140個甲基化位點納入本次研究。對mRNA數(shù)據(jù)的表達(dá)值進(jìn)行l(wèi)og2(count+1)轉(zhuǎn)化，由于不同組學(xué)存在異質(zhì)性問題，將轉(zhuǎn)化后的mRNA數(shù)據(jù)以及甲基化數(shù)據(jù)進(jìn)行Z-score標(biāo)準(zhǔn)化處理。

表1 基于NMI評價指標(biāo)對五種無監(jiān)督聚類方法結(jié)果比較

對預(yù)處理得到的數(shù)據(jù)分別使用ANF、SNF、K-means、iCluster以及系統(tǒng)聚類方法分析，提示聚類分為兩類，因此將卵巢癌鉑類藥物敏感患者分為兩個免疫亞型即免疫I型和免疫II型，通過log-rank檢驗比較亞型間的生存差異。為增加方法間的可比性，對于SNF和ANF兩種方法，分別選擇log-rank檢驗得到P值最小時的參數(shù)k，其中ANF中k分別選擇30和15。結(jié)果參見表2。

表2 無監(jiān)督聚類得到免疫分子分型的log-rank生存分析結(jié)果比較

結(jié)果顯示對ANF得到的無監(jiān)督聚類標(biāo)簽進(jìn)行l(wèi)og-rank檢驗P<0.001(HR=2.900，95%CI(1.809，4.648))，P值小于其他四個方法的P值，其中通過系統(tǒng)聚類進(jìn)行生存分析模型不收斂，同時多組學(xué)整合無監(jiān)督聚類結(jié)果優(yōu)于單組學(xué)無監(jiān)督聚類結(jié)果。圖1表示利用以上方法進(jìn)行多組學(xué)整合得到的Kaplan-Meier生存曲線。

圖1-A顯示利用ANF方法得到的免疫分型I型和II型的生存曲線，發(fā)現(xiàn)I型預(yù)后好于II型預(yù)后。為了進(jìn)一步證明多組學(xué)聚類效果優(yōu)于單組學(xué)聚類效果，使用ANF方法得到的mRNA、甲基化以及mRNA和甲基化多組學(xué)整合的聚類結(jié)果計算輪廓系數(shù)，參見圖2，A～C圖中縱軸表示每個樣本及其所在的亞組，橫軸表示每個樣本的輪廓系數(shù)值，最終取平均得到該組學(xué)的輪廓系數(shù)值，值越接近1表示聚類效果越合理。結(jié)果顯示，通過多組學(xué)聚類的輪廓系數(shù)0.42高于單組學(xué)聚類輪廓系數(shù)(mRNA：0.32；甲基化：0.3)，進(jìn)一步說明多組學(xué)在無監(jiān)督聚類亞型分析中的價值。

圖1 多組學(xué)整合的化療敏感患者免疫分型生存曲線

為了控制其他因素對生存的影響，將年齡(≥60歲和<60歲)，化療緩解情況(完全緩解和非完全緩解)，腫瘤分化程度(高分化G1，中分化G2，低分化G3，未分化G4)，腫瘤臨床分期(I期，II期，III期，IV期)列為協(xié)變量，Cox比例風(fēng)險回歸模型結(jié)果見表3。

結(jié)果表明，校正上述協(xié)變量后，I型預(yù)后仍好于II型預(yù)后。可以看出ANF得到的免疫亞型是卵巢癌化療敏感人群的獨立預(yù)后因素。

為了進(jìn)一步確定I型和II型的差異變量，本研究使用LASSO算法進(jìn)行組間差異變量篩選，得到43個差異mRNA以及39個差異甲基化位點，將得到的差異變量在不同樣本間表達(dá)量尺度歸一化，使用R包pheatmap畫出熱圖，表示如圖3。

表3 調(diào)整協(xié)變量后的Cox比例風(fēng)險回歸模型結(jié)果

根據(jù)熱圖可以看出不同亞型的組間差別。據(jù)文獻(xiàn)報道，BMPR1B在乳腺癌中高表達(dá)與較差的生存預(yù)后相關(guān)[8]。PIK3CD在其他癌癥如結(jié)腸癌組織和CRC細(xì)胞系中過表達(dá)，是結(jié)腸癌患者總體生存的獨立預(yù)測因子，并且PIK3CD與CILA4和原發(fā)性免疫缺陷有關(guān)，有理由認(rèn)為BMPR1B和PIK3CD的高表達(dá)可能與卵巢癌化療敏感的患者生存預(yù)后較差相關(guān)[9-10]。甲基化位點中cg16097079，cg17475918分別對應(yīng)HLA-C、HLA-E，屬于主要組織相容性復(fù)合體中(MHC)的I類分子，其表達(dá)的喪失提供了關(guān)鍵的免疫逃避策略[11]。巨噬細(xì)胞遷移抑制因子(MIF)是一種炎性介質(zhì)，可促進(jìn)炎癥細(xì)胞因子的產(chǎn)生，并具有趨化特性，可招募巨噬細(xì)胞和T細(xì)胞，目前有研究探討可將MIF作為靶標(biāo)對卵巢癌以及泌尿生殖系統(tǒng)癌癥患者進(jìn)行治療，本研究中在甲基化位點cg20377673所對應(yīng)的基因以及mRNA中均發(fā)現(xiàn)MIF表達(dá)在得到的免疫亞型間存在差異[12-14]。因此，針對鉑類藥物敏感的卵巢癌患者的不同免疫亞型，使用特定的免疫靶點治療可能增加患者的預(yù)后生存。

圖3 不同分子亞型中差異基因表達(dá)的熱圖

討 論

與傳統(tǒng)的使用單組學(xué)數(shù)據(jù)聚類發(fā)現(xiàn)疾病亞型相比，包含多種組學(xué)信息可以彌補(bǔ)單組學(xué)的局限性，生物學(xué)解釋更合理，有效提高患者預(yù)后效果。利用樣本多組學(xué)數(shù)據(jù)進(jìn)行網(wǎng)絡(luò)構(gòu)建的ANF算法在原始SNF算法的基礎(chǔ)上簡化了迭代過程，使計算上更為簡便且效果良好，提高運算效率。通過矩陣變換使其在解釋上更為合理，同時可對不同組學(xué)信息進(jìn)行加權(quán)，更符合生物學(xué)意義。

模擬實驗結(jié)果提示通過改變差異變量的均數(shù)差值以及差異變量百分比，當(dāng)亞型總體差異較小時，ANF方法聚類效果明顯優(yōu)于其他方法。均數(shù)差異越大，差異變量百分比越高時，NMI越高，聚類效果越好，當(dāng)不同亞型總體差異足夠大時，除系統(tǒng)聚類外其他幾種方法效果趨于一致。

由于鉑類藥物敏感的卵巢癌患者總生存較低，有必要對這類患者進(jìn)行進(jìn)一步治療提高患者生存。本次研究發(fā)現(xiàn)對鉑類藥物敏感的卵巢癌患者的新的免疫亞型，且該免疫亞型中通過LASSO篩選得到的差異mRNA和差異甲基化位點中，有目前報道的免疫靶點MIF等，因此針對不同免疫亞型，選擇免疫靶點進(jìn)行治療可能是鉑類藥物化療后的又一治療策略。

由于收集多組學(xué)數(shù)據(jù)較為困難，目前大多數(shù)利用多組學(xué)得到的分子亞型在驗證上難度較大。對于無監(jiān)督聚類，生存分析的結(jié)局指標(biāo)常用來對獲得亞型進(jìn)行評價分析。方法中參數(shù)的選擇可能會對假設(shè)檢驗中的P值產(chǎn)生較大影響，因此這并不能作為唯一的評價指標(biāo)，在未知組學(xué)權(quán)重的情況下，權(quán)重的設(shè)置通常較難選擇，實際中可以根據(jù)已有先驗知識進(jìn)行設(shè)置。