寇傳華 張媛媛

（青島理工大學(xué)信息與控制工程學(xué)院 青島 266520）

1 引言

人類復(fù)雜疾病的發(fā)生和發(fā)展通常不是由單一因素引起的，而是由多種因素相互作用引起的功能性障礙［1~2］。因此，基于網(wǎng)絡(luò)模型挖掘復(fù)雜疾病生物標(biāo)志物是復(fù)雜疾病相關(guān)研究的重要內(nèi)容［3~4］。然而，目前基于網(wǎng)絡(luò)模型挖掘疾病標(biāo)記物的研究中很少涉及到對所構(gòu)建網(wǎng)絡(luò)模型真實(shí)性的評估。

DNA 甲基化作為一種重要的表觀遺傳修飾，在生物體的生命過程中，環(huán)境會影響DNA 甲基化水平從而調(diào)控基因表達(dá)出不同的表型。因此，基于甲基化數(shù)據(jù)識別癌癥致病因子有助于在早期發(fā)現(xiàn)疾病并對疾病治療方案給出建議［5］。同時(shí)，研究基因間甲基化水平變化的相關(guān)性可以在癌癥早期階段找到與癌癥相關(guān)的基因模塊。目前，基于DNA甲基化數(shù)據(jù)并利用其數(shù)據(jù)特點(diǎn)，出現(xiàn)了大量基于網(wǎng)絡(luò)模型的直接關(guān)聯(lián)所有CpG 位點(diǎn)對的方法［6~7］。因?yàn)槊總€(gè)基因?qū)?yīng)多個(gè)甲基化位點(diǎn)并且不同基因?qū)?yīng)的甲基化位點(diǎn)數(shù)目是不同的，所以構(gòu)建基因網(wǎng)絡(luò)的關(guān)鍵在于如何把CpG 位點(diǎn)的甲基化水平映射到基因?qū)用嫔?。根?jù)不同的融合方法，可以將其分為基于平均值的方法［6］、總體值的方法［7］和差異得分的方法［8］。Ma 等［6］在DNA 甲基化網(wǎng)絡(luò)的構(gòu)建中，將基因中一些位點(diǎn)的平均甲基化水平作為基因的甲基化水平，用Pearson 相關(guān)系數(shù)（PCC）來測量基因之間的相關(guān)性。Bartlett 等［7］基于典型相關(guān)分析（CCA）的方法構(gòu)建甲基化基因網(wǎng)絡(luò)來識別與癌癥相關(guān)的模塊；West 等［8］提出的EpiMod 方法，該方法將表型上每對基因差異得分（diffscore）的平均值作為邊緣權(quán)重來識別與年齡相關(guān)的甲基化模塊。綜上所述，根據(jù)DNA 甲基化數(shù)據(jù)的特點(diǎn)，在基因?qū)用嫔喜捎貌煌瑴y度的多位點(diǎn)甲基化融合方法可以得到不同的基因網(wǎng)絡(luò)［9~10］，不同網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)會影響疾病相關(guān)生物標(biāo)志物的識別［11］。

在本文中，受Nguyen 等提出擾動聚類算法的啟發(fā)［12］，我們提出了ENCPer 方法。該方法評估用不同的方法構(gòu)建的基于DNA 甲基化數(shù)據(jù)的基因網(wǎng)絡(luò)的穩(wěn)定性（圖1）。我們實(shí)驗(yàn)的基本過程如下：首先，基于真實(shí)的DNA 甲基化數(shù)據(jù)用不同的方法構(gòu)建基因網(wǎng)絡(luò)。其次，將高斯噪聲數(shù)據(jù)作為擾動加入到真實(shí)的DNA 甲基化數(shù)據(jù)中，最后，重建基于擾動數(shù)據(jù)的基因網(wǎng)絡(luò)。為衡量不同網(wǎng)絡(luò)構(gòu)建方法的優(yōu)劣，本文提出了一種網(wǎng)絡(luò)穩(wěn)定性度量方法ENCPer，就如何基于DNA甲基化數(shù)據(jù)構(gòu)建網(wǎng)絡(luò)提供建議。

圖1 本文框架

2 材料與方法

2.1 基因加權(quán)網(wǎng)絡(luò)的構(gòu)建

在獲得上述差異甲基化基因后，本文從三個(gè)角度：融合數(shù)據(jù)的平均值，總體值和差異得分構(gòu)建了基因網(wǎng)絡(luò)。在本節(jié)中描述了三種融合角度下基于PCC，CCA和diffscore基因網(wǎng)絡(luò)的構(gòu)建方法。

基于PCC方法的基因加權(quán)網(wǎng)絡(luò)的構(gòu)建：對于每個(gè)基因，將一個(gè)基因中所有甲基化位點(diǎn)的平均甲基化水平作為該基因的甲基化水平。對于正常樣本和癌癥樣本，使用PCC 計(jì)算兩個(gè)基因之間的相關(guān)性，并分別獲得處于正常和癌癥狀態(tài)的兩個(gè)基因網(wǎng)絡(luò)鄰接矩陣和，然后使用正常和癌癥網(wǎng)絡(luò)對應(yīng)邊之差得到差分網(wǎng)絡(luò)鄰接矩陣ΔMpcc。

基于CCA 方法的基因加權(quán)網(wǎng)絡(luò)的構(gòu)建：基因X和Y包含不同數(shù)量的甲基化位點(diǎn)，用線性組合U1=aT?X和V1=bT?Y描述基因。兩個(gè)基因的相關(guān)性由下面公式描述：

cov(U1,V1)是U1和V1的協(xié)方差，var(U1) 和var(V1) 分別是U1和V1的方差。在所有可能的線性組合中找出使ρ最大的線性組合的對U1和V1，通過相關(guān)性ρ建立基因網(wǎng)絡(luò)，這三個(gè)網(wǎng)絡(luò)的鄰接矩陣分別表示為，和ΔMcca。

基于diffscore 方法的基因加權(quán)網(wǎng)絡(luò)的構(gòu)建：可以使用某種距離度量（歐氏距離或相對熵）比較腫瘤樣本與正常樣本之間每個(gè)位點(diǎn)的差異，并選擇中位數(shù)或者平均值作為基因的差異得分（這里使用平均值和中位數(shù)），表示為GDS(X)。通過計(jì)算基因X和Y之間的權(quán)重wXY構(gòu)造diffscore 差分網(wǎng)絡(luò)，diffscore 網(wǎng)絡(luò)的鄰接矩陣表示為ΔMdiffscore。公式定義如下：

2.2 顯著性計(jì)算

本文使用相似性置換的方法來顯示差分網(wǎng)絡(luò)中加權(quán)邊緣的重要性。首先，對于原始數(shù)據(jù)中的樣本標(biāo)簽進(jìn)行nperm次置換。其次，每次置換后計(jì)算正常樣本和癌癥樣本對應(yīng)的基因網(wǎng)絡(luò)和差分網(wǎng)絡(luò)，差分網(wǎng)絡(luò)鄰接矩陣表示為perm_ΔMmethod。最后，使用下面的公式計(jì)算每個(gè)邊緣權(quán)重的p值。

2.3 不同網(wǎng)絡(luò)的穩(wěn)定性評價(jià)

本文提出了ENCPer 算法，在原始的甲基化數(shù)據(jù)中隨機(jī)添加高斯噪聲θ，并利用上述三種計(jì)算基因間權(quán)重的方法重新構(gòu)造基于擾動數(shù)據(jù)構(gòu)建的擾動基因網(wǎng)絡(luò)。并通過以下定義的累積分布函數(shù)（CDF）來評估由這三種方法構(gòu)建的網(wǎng)絡(luò)。

公式中I()?是一個(gè)指標(biāo)函數(shù)，E是一組邊，wij代表基因i和基因j在不同方法構(gòu)建的真實(shí)網(wǎng)絡(luò)中網(wǎng)絡(luò)邊緣的權(quán)重；同樣地，代表基因i和基因j在不同方法構(gòu)建的擾動網(wǎng)絡(luò)中網(wǎng)絡(luò)邊緣的權(quán)重。τ是CDF 的閾值，num(E)代表相應(yīng)網(wǎng)絡(luò)中的邊數(shù)。具有最高穩(wěn)定性的網(wǎng)絡(luò)將被認(rèn)為是最接近真實(shí)生物系統(tǒng)的網(wǎng)絡(luò)。

3 結(jié)果

3.1 基于不同方法的網(wǎng)絡(luò)特征分析

使用PCC［6］，CCA［7］和diffscore［8］方法構(gòu)建在BRCA和COAD數(shù)據(jù)集上的基因加權(quán)網(wǎng)絡(luò)［11］。根據(jù)權(quán)度（degree），緊密度（close）和page rank 可以比較在不同網(wǎng)絡(luò)中對應(yīng)節(jié)點(diǎn)的重要性和不同網(wǎng)絡(luò)的不同特征。通過對不同網(wǎng)絡(luò)中的邊重疊數(shù)目的比較，可以看出根據(jù)不同方法構(gòu)建的基因網(wǎng)絡(luò)差異性很大（圖2）。根據(jù)不同網(wǎng)絡(luò)中節(jié)點(diǎn)重要性得分的排序選擇前500個(gè)基因進(jìn)行重疊比較（圖3）。發(fā)現(xiàn)除了在BRCA 數(shù)據(jù)集上對節(jié)點(diǎn)進(jìn)行基于權(quán)度評估中三個(gè)網(wǎng)絡(luò)的交集較大，其他情況下網(wǎng)絡(luò)的交集較少。這表明用不同方法構(gòu)建網(wǎng)絡(luò)對結(jié)果的影響很大，所以選擇合適的網(wǎng)絡(luò)構(gòu)建方法對于了解真實(shí)的生物系統(tǒng)或識別與疾病真正相關(guān)的標(biāo)記具有重要意義。

圖2 三種不同網(wǎng)絡(luò)中的邊的交集

圖3 在不同節(jié)點(diǎn)重要性評估中三個(gè)網(wǎng)絡(luò)的重疊

3.2 網(wǎng)絡(luò)穩(wěn)定性分析

在BRCA 和CODA 數(shù)據(jù)集上，分別基于PCC，CCA 和diffscore 方法構(gòu)建真實(shí)基因網(wǎng)絡(luò)。在這里，我們選擇兩個(gè)距離指標(biāo)，歐氏距離（Eu）［13］和相對熵（RE）［14］以及均值（mean），中位數(shù)（median）兩個(gè)綜合指標(biāo)來計(jì)算基因的差異得分。我們將這四種條件（Eu_mean，Eu_median，RE_mean，RE_median）應(yīng)用于diffscore 并于PCC 和CCA 方法構(gòu)建的基因網(wǎng)絡(luò)進(jìn)行比較，我們向?qū)嶋H數(shù)據(jù)中添加了不同的噪聲θ=(0.01 ,0.05,0.1,0.2,0.3,0.5 )。根據(jù)式（4）對CDF進(jìn)行計(jì)算，通過改變閾值可以得到不同噪聲下的四個(gè)基因網(wǎng)絡(luò)的穩(wěn)定性得分（圖4）。

圖4 BRCA數(shù)據(jù)集（左側(cè)）和COAD數(shù)據(jù)集（右側(cè)）上的三個(gè)網(wǎng)絡(luò)在噪聲和閾值下的穩(wěn)定性評估

通過觀察不同網(wǎng)絡(luò)在不同數(shù)據(jù)集中的性能，我們發(fā)現(xiàn)當(dāng)τ=0 時(shí)，基于CCA 方法構(gòu)建的基因網(wǎng)絡(luò)穩(wěn)定性最差。這可能與以下原因有關(guān)：CCA的權(quán)重計(jì)算方法綜合考慮了所有位點(diǎn)的信息，錯(cuò)誤地合成了真實(shí)信號和噪聲，因此對隨機(jī)添加的噪聲敏感。一方面，基于diffscore 方法構(gòu)建的網(wǎng)絡(luò)無論使用哪種距離測度和綜合指標(biāo)都具有中等穩(wěn)定性?；赑CC 的網(wǎng)絡(luò)雖然在COAD 數(shù)據(jù)集中具有較高的穩(wěn)定性，但在BRCA 數(shù)據(jù)集中穩(wěn)定性較低。另一方面，在不同的噪聲下，基于PCC 和CCA 方法構(gòu)建的網(wǎng)絡(luò)具有魯棒性，基于diffscore 的網(wǎng)絡(luò)除了根據(jù)距離測度和綜合指標(biāo)選擇出來的相對熵和均值具有魯棒性外，其他不具有魯棒性。因此，結(jié)合以上兩個(gè)方面的分析，基于PCC 的網(wǎng)絡(luò)和基于diffscore 的網(wǎng)絡(luò)具有更高的穩(wěn)定性，并且對噪聲的魯棒性更高。

3.3 疾病相關(guān)生物標(biāo)志物的鑒定

在得到BRCA 和COAD 數(shù)據(jù)集中表現(xiàn)最穩(wěn)定的網(wǎng)絡(luò)后，使用R 語言igraph 包中的cluster_fast_greedy 算法優(yōu)化模塊化得分來得到稠密子圖。使用clusterProfiler 包對兩個(gè)網(wǎng)絡(luò)中的子圖進(jìn)行KEGG 和GO 富集分析。結(jié)果表明，子圖中的基因在豐富細(xì)胞分化，細(xì)胞命運(yùn)等與疾病相關(guān)的GO 術(shù)語的同時(shí)也豐富了多種疾病的相關(guān)通路，例如細(xì)胞粘附分子，胃癌和乳腺癌作用等。為了展顯最穩(wěn)定網(wǎng)絡(luò)的性能，我們與其他網(wǎng)絡(luò)富集分析的結(jié)果進(jìn)行比較，實(shí)驗(yàn)結(jié)果表明，在COAD 數(shù)據(jù)集上基于diff?score 方法構(gòu)建的基因網(wǎng)絡(luò)中找不到明顯的稠密子圖；在基于PCC方法的基因網(wǎng)絡(luò)中發(fā)現(xiàn)了疾病相關(guān)富集通路（如皮膚癌、乳腺癌等），但這些富集通路在基于CCA方法的網(wǎng)絡(luò)中沒有發(fā)現(xiàn)。

4 結(jié)語

DNA 甲基化是基因組與環(huán)境之間的聯(lián)系，通常在癌癥發(fā)生的早期階段發(fā)生變化［15］。精確鑒定出與甲基化變化相關(guān)的生物標(biāo)記物可作為早期癌癥風(fēng)險(xiǎn)評估的風(fēng)險(xiǎn)因素。此外，由于癌癥是多種因素相互作用的結(jié)果，所以，構(gòu)建基于DNA 甲基化的基因網(wǎng)絡(luò)對于識別癌癥相關(guān)功能模塊具有重要意義。在兩個(gè)真實(shí)數(shù)據(jù)集BRCA 和COAD 中，我們分別從平均值，總體值和差異得分三種角度使用三種不同的方法構(gòu)建基因加權(quán)網(wǎng)絡(luò)。使用ENCPer方法評估構(gòu)建出不同網(wǎng)絡(luò)的穩(wěn)定性。通過比較這些網(wǎng)絡(luò)特性，我們發(fā)現(xiàn)不同方法構(gòu)建的網(wǎng)絡(luò)有很大的差異，這也充分說明了研究網(wǎng)絡(luò)穩(wěn)定性評估策略對于后續(xù)工作的重要性。通過對不同網(wǎng)絡(luò)的穩(wěn)定性評估，我們發(fā)現(xiàn)在BRCA數(shù)據(jù)集中，基于diffscore方法構(gòu)建的網(wǎng)絡(luò)最穩(wěn)定，而在COAD 數(shù)據(jù)集中，基于PCC 方法構(gòu)建的網(wǎng)絡(luò)最穩(wěn)定。針對這兩個(gè)穩(wěn)定的網(wǎng)絡(luò)，通過社區(qū)檢測算法分別獲得了兩個(gè)網(wǎng)絡(luò)的稠密子圖，并對獲得的稠密子圖進(jìn)行了GO 和KEGG富集分析，得到的富集分析結(jié)果顯示了網(wǎng)絡(luò)結(jié)構(gòu)的生物學(xué)意義以及基于擾動算法構(gòu)建出穩(wěn)定網(wǎng)絡(luò)的有效性。