基于DNA 甲基化數(shù)據(jù)的擾動算法評估網(wǎng)絡構建穩(wěn)定性?

寇傳華 張媛媛

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

1 引言

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

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

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

圖1 本文框架

2 材料與方法

2.1 基因加權網(wǎng)絡的構建

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

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

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

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

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

2.2 顯著性計算

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

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

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

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

3 結果

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

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

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

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

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

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

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

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

3.3 疾病相關生物標志物的鑒定

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

4 結語

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