馬欣宇，王 嬋，胡躍清

(復(fù)旦大學(xué) 生命科學(xué)學(xué)院 生物統(tǒng)計(jì)學(xué)與計(jì)算生物學(xué)系，上海 200433)

基于病例對照組的檢驗(yàn)差異甲基化基因功能區(qū)域的方法

馬欣宇，王 嬋，胡躍清

(復(fù)旦大學(xué) 生命科學(xué)學(xué)院 生物統(tǒng)計(jì)學(xué)與計(jì)算生物學(xué)系，上海 200433)

DNA甲基化是一種重要的表觀遺傳學(xué)修飾，能夠參與基因表達(dá)調(diào)控從而影響生理活動(dòng)．在病例-對照組研究中，利用甲基化水平定位與疾病相關(guān)聯(lián)的基因功能區(qū)域，有助于了解甲基化調(diào)控基因表達(dá)的機(jī)制，為后續(xù)研究提供疾病關(guān)聯(lián)基因的線索．目前檢驗(yàn)甲基化差異區(qū)域的方法存在一定不足： 一是缺乏定義差異甲基化區(qū)域的統(tǒng)一標(biāo)準(zhǔn)；二是檢測結(jié)果存在橫跨多個(gè)基因功能區(qū)域的現(xiàn)象，難以進(jìn)行生物學(xué)解釋．本文基于得分檢驗(yàn)(Z-score test)提出了以基因功能區(qū)域?yàn)閱挝坏臋z驗(yàn)差異甲基化水平的方法cZST．該方法的特點(diǎn)是依靠位點(diǎn)物理距離和染色體全局信息修正協(xié)方差矩陣．模擬研究表明cZST在亞硫酸氫鹽測序數(shù)據(jù)樣本量小、讀數(shù)低、擾動(dòng)位點(diǎn)多的情形下具有高功效，且優(yōu)于現(xiàn)有方法．將cZST運(yùn)用到乳腺導(dǎo)管癌真實(shí)數(shù)據(jù)中，檢測出與疾病顯著相關(guān)的661個(gè)基因功能區(qū)和81個(gè)單位點(diǎn)．基因注釋分析說明了cZST方法的有效性，并發(fā)現(xiàn)Dbl、Pleckstrin基因家族及Rho調(diào)控基因的甲基化與乳腺導(dǎo)管癌緊密相關(guān)，為后續(xù)研究提供指引．

差異甲基化檢驗(yàn)； 亞硫酸氫鹽測序； 功能富集分析

DNA甲基化是發(fā)生在生物DNA鏈核苷酸上的一種化學(xué)修飾．在真核細(xì)胞中，S-腺苷甲硫氨酸上的甲基會在DNA甲基轉(zhuǎn)移酶的作用下，轉(zhuǎn)移到DNA分子CpG序列的胞嘧啶上[1]．人類DNA中約70%～90%CpG二核苷酸存在甲基化修飾[2]，并且甲基化的CpG二核苷酸大部分聚集在基因的5’端，稱為CpG島，長度在幾百至一千個(gè)堿基對(bp)左右，它們大多分布在基因的上游、啟動(dòng)子和外顯子區(qū)域中[3]．

研究表明CpG島上的甲基化會顯著影響基因的轉(zhuǎn)錄調(diào)節(jié)，通常會使轉(zhuǎn)錄失活．例如，在人和嚙齒類動(dòng)物中，病毒序列可以與插入的基因交聯(lián)而甲基化，從而使該基因沉默[4]．DNA甲基化修飾是可逆的，并且不直接依賴于序列，能夠參與基因表達(dá)調(diào)控[5-6]．因此，DNA甲基化作為研究人類疾病、發(fā)育、衰老的重要媒介，已成為當(dāng)今基因組學(xué)的重要研究對象之一．

第二代測序技術(shù)，由于其高靈敏度和特異性，已經(jīng)成為表觀基因組研究的有力工具．BS-seq(Bisulfite Sequencing)使用亞硫酸氫鹽處理DNA，未甲基化的胞嘧啶會轉(zhuǎn)化為尿嘧啶，而甲基化的胞嘧啶不發(fā)生變化[7]．因此，BS-seq能獲得基因組的全堿基甲基化信息，從而計(jì)算出單核苷酸分辨率下甲基化在細(xì)胞樣本中的發(fā)生率，即甲基化水平值，它等于被甲基化信號強(qiáng)度與總信號強(qiáng)度之比．由于全基因組BS-seq較為昂貴，許多研究者使用簡易BS-seq(RRBS)，它基因組覆蓋率不到10%，測序成本相對較低[8-9]．

隨著芯片技術(shù)和高通量測序技術(shù)的發(fā)展，全表觀基因組關(guān)聯(lián)研究成為一個(gè)研究熱點(diǎn)．它通過檢驗(yàn)基因組上，尤其是CpG島上，DNA核苷酸甲基化的分布差異，搜索與性狀關(guān)聯(lián)性強(qiáng)的位點(diǎn)或區(qū)域．已有的方法包括COHCAP[10]、BISEQ[11]、HMM-Fisher[12]等，它們能夠較為準(zhǔn)確地檢測BS-seq單位點(diǎn)的甲基化差異，但存在一定的局限性．由于多數(shù)CpG位點(diǎn)成簇存在，共同影響所在基因，而一個(gè)基因功能區(qū)域內(nèi)各甲基化位點(diǎn)的表達(dá)調(diào)控作用是互補(bǔ)累加的，因此單位點(diǎn)差異甲基化檢驗(yàn)提供的信息往往是不充分甚至片面的．而針對連續(xù)多個(gè)CpG位點(diǎn)的差異甲基化區(qū)域檢驗(yàn)，缺乏劃分或拼接區(qū)域的統(tǒng)一標(biāo)準(zhǔn)[13]，難以比較和重復(fù)驗(yàn)證．同時(shí)其分析結(jié)果常常橫跨多個(gè)功能區(qū)域甚至兩個(gè)基因，這給進(jìn)行生物學(xué)解釋帶來了困難．

本文提出一種基于得分的差異甲基化檢驗(yàn)，它考慮了差異的方向性，量化了基因功能區(qū)域的差異顯著程度，其生物學(xué)意義明顯，而且方便后續(xù)研究．同時(shí)針對甲基化BS-seq數(shù)據(jù)樣本量較小、測序讀數(shù)低的特點(diǎn)，提出了依靠位點(diǎn)物理距離和染色體全局信息修正協(xié)方差矩陣的策略，它顯著地提高了檢驗(yàn)功效(真陽性率,TPR)．

1 符號與方法

(1)

在原假設(shè)H0下，即病例組和對照組的甲基化水平一致時(shí)，U近似服從均值為0，協(xié)方差矩陣為V的多維正態(tài)分布，從而構(gòu)造檢驗(yàn)統(tǒng)計(jì)量S=UTV-1U，其中

(2)

由于測序成本高，BS-seq數(shù)據(jù)樣本量小，測序讀數(shù)低．而低讀數(shù)位點(diǎn)會導(dǎo)致甲基化水平的估計(jì)值誤差較大，從而協(xié)方差偏差也較大，進(jìn)而影響檢驗(yàn)功效．為了彌補(bǔ)這種不足，我們引入全局甲基化信息，利用同一染色體上其他位點(diǎn)甲基化水平的相關(guān)系數(shù)，來校正低讀數(shù)位點(diǎn)協(xié)方差的估計(jì)，從而增加檢驗(yàn)功效．

具體校正策略如下．先計(jì)算整條染色體上兩兩位點(diǎn)的物理間距以及相應(yīng)的甲基化水平的相關(guān)系數(shù)(Correlation coefficient)，并進(jìn)行局部加權(quán)線性回歸獲得LOWESS曲線，作為全局信息以備后面使用．然后設(shè)定兩個(gè)讀數(shù)閾值T1=12和T2=4，對位點(diǎn)進(jìn)行分類．給定某位點(diǎn)，記T為所有個(gè)體的總信號讀數(shù)(甲基化和非甲基化信號之和)的平均值．若T≤T2，則認(rèn)為此位點(diǎn)的數(shù)據(jù)質(zhì)量低，將其從數(shù)據(jù)集中剔除；如果T2

該修正策略基于甲基化BS-seq數(shù)據(jù)如下的特征： 同一條染色體上相同物理距離的兩位點(diǎn)甲基化水平的相關(guān)系數(shù)的波動(dòng)程度小，因此這些相關(guān)系數(shù)的平均值與低讀數(shù)位點(diǎn)真實(shí)相關(guān)系數(shù)較為接近．我們選用一個(gè)RRBS真實(shí)數(shù)據(jù)集[15]，對1號染色體上位點(diǎn)間甲基化水平相關(guān)系數(shù)進(jìn)行統(tǒng)計(jì)分析，其平均值及變異系數(shù)(Coefficient of Variation, CV)相對物理距離的散點(diǎn)圖見圖1．我們發(fā)現(xiàn)在點(diǎn)間距小于一定范圍(如400bp)的兩位點(diǎn)對應(yīng)β值相關(guān)系數(shù)明顯大于0且與間距呈反比，并且當(dāng)間距近時(shí)，同一物理距離的兩位點(diǎn)β值相關(guān)系數(shù)的變異系數(shù)很低，說明相關(guān)系數(shù)波動(dòng)程度小，這說明了修正策略的合理性．

此外，如果忽略任意兩位點(diǎn)之間的甲基化水平的相關(guān)性，即相當(dāng)于令V的非對角元素全部為零，對角元為樣本方差，這種得分檢驗(yàn)記作nZST(non-related ZST)，將它與cZST、ZST一起進(jìn)行比較．

統(tǒng)計(jì)量S的p值可采取置換[16]來獲?。诿看沃脫Q時(shí)，我們隨機(jī)打亂每個(gè)個(gè)體的患病狀態(tài)，然后計(jì)算新統(tǒng)計(jì)量Sm，1≤m≤M．則p值為

(3)

在模擬研究中，基于R次重復(fù)模擬數(shù)據(jù)的p值pr，r=1,2,…,R，在顯著性水平α下，該方法的功效RTP為

(4)

圖1 β值相關(guān)系數(shù)的均值和變異系數(shù)散點(diǎn)圖Fig.1 Mean value and CV of β values’ correlation coefficient(a) RRBS數(shù)據(jù)1號染色體上同距離兩位點(diǎn)的β值相關(guān)系數(shù)的均值隨物理間距的變化情況；(b) 相關(guān)系數(shù)的變異系數(shù)隨物理間距的變化情況．實(shí)線為LOWESS擬合平滑曲線．

2 隨機(jī)模擬研究

為了模擬實(shí)際BS-seq數(shù)據(jù)特征，我們通過多維正態(tài)分布刻畫位點(diǎn)間甲基化水平相關(guān)性，通過二項(xiàng)分布模擬產(chǎn)生測序信號．假設(shè)某區(qū)域的總長度為800bp，含L個(gè)CpG位點(diǎn)，隨機(jī)從1到800中選取L個(gè)不同的位置記為D1,D2,…,DL．考慮兩個(gè)多維正態(tài)分布N(μ1,Σ)和N(μ0,Σ)，μ1、μ0為兩個(gè)均值向量，對應(yīng)協(xié)方差矩陣Σ=(σil)L×L，其中的對角元均為1，非對角元為

σil=exp(-dil/d0),

(5)

其中dil=|Di-Dl|，d0為默認(rèn)任意兩位點(diǎn)β值存在相關(guān)性的距離參數(shù)，在我們的隨機(jī)模擬中取d0=400bp．該隨機(jī)模型建立了物理距離與位點(diǎn)相關(guān)性之間的聯(lián)系，距離更近的兩個(gè)位點(diǎn)甲基化水平擁有更大的相關(guān)系數(shù)，反之亦然．

以多維正態(tài)分布N(μ1,Σ)或N(μ0,Σ)產(chǎn)生隨機(jī)數(shù)，通過逆Logistic變換得到取值在[0,1]的L個(gè)值P1,P2,…,PL，作為觀察到甲基化信號的理論概率；設(shè)單位點(diǎn)測序總讀數(shù)的期望均為T0，通過泊松分布P(T0)產(chǎn)生L個(gè)隨機(jī)數(shù)N1,N2,…,NL,作為測序總讀數(shù)．最終通過二項(xiàng)分布B(Nl,Pl)產(chǎn)生L個(gè)位點(diǎn)上的甲基化和未甲基化讀數(shù)，進(jìn)而得到甲基化水平值β1,β2,…,βL．因此，每次模擬產(chǎn)生L個(gè)不同的位置，從N(μ1,Σ)出發(fā)可以產(chǎn)生n1個(gè)病例的L個(gè)甲基化水平值，從N(μ0,Σ)出發(fā)可產(chǎn)生n0個(gè)對照的L個(gè)甲基化水平值．取n1=n0=6．

實(shí)驗(yàn)1： 在測序讀數(shù)足夠大且位點(diǎn)存在相關(guān)性的條件下，比較ZST、cZST和nZST 3種檢驗(yàn)方法的表現(xiàn)．模擬研究中cZST方法全部使用修正的協(xié)方差．基于原假設(shè)，取μ1=μ0=0，產(chǎn)生4000組無差異數(shù)據(jù)；基于備擇假設(shè)，取μ1≠μ0，產(chǎn)生1000組有差異數(shù)據(jù)．假設(shè)這5000組數(shù)據(jù)對應(yīng)的區(qū)域位于同一染色體上，且兩兩區(qū)域間距較遠(yuǎn)．分別使用上述3種方法進(jìn)行檢驗(yàn)，將P值從小到大排序，依次選取1到前2000個(gè)作為顯著性區(qū)域，計(jì)算錯(cuò)誤發(fā)現(xiàn)率(False Discovery Rate, FDR)．

當(dāng)T0=25，即測序讀數(shù)較大時(shí)；或令β值完全等于甲基化理論概率，即βi=Pi，i=1,2,…,L時(shí)，分別產(chǎn)生模擬數(shù)據(jù)，F(xiàn)DR結(jié)果如圖2所示．圖2結(jié)果表明，利用全局甲基化水平的相關(guān)系數(shù)信息的cZST控制FDR的能力遠(yuǎn)強(qiáng)于使用原始協(xié)方差矩陣的ZST及不考慮位點(diǎn)相關(guān)性的nZST．

圖2 3種方法在測序總讀數(shù)大時(shí)FDR比較Fig.2 FDR values of three methods given a large reading depth3種方法選取不同顯著位點(diǎn)個(gè)數(shù)時(shí)對應(yīng)的FDR值．μ0始終為0，原假設(shè)下μ1=0，備擇假設(shè)下μ1=(0,0,0,1,1,1,0,0)．

實(shí)驗(yàn)2： 當(dāng)位點(diǎn)之間的β值相關(guān)性較弱時(shí)，全局甲基化相關(guān)系數(shù)提供的真實(shí)信息減少，cZST可能受到干擾，探究該情況下cZST的穩(wěn)健性．在模型中加入一個(gè)用于調(diào)整相關(guān)性的權(quán)重系數(shù)c，將多維正態(tài)分布協(xié)方差矩陣設(shè)為

Σ′=cΣ+(1-c)I,c∈[0,1]．

則c取值越小，位點(diǎn)間模擬β值的相關(guān)性越弱．其余參數(shù)不變，原假設(shè)下產(chǎn)生數(shù)據(jù)計(jì)算3種方法的第一類錯(cuò)誤率(type I error rate)，備擇假設(shè)下產(chǎn)生數(shù)據(jù)比較它們功效的差異，α=0.05，結(jié)果如圖3所示．在三者均能控制住第一類錯(cuò)誤率的前提下，c>0.5時(shí)，cZST功效優(yōu)勢持續(xù)增加；cZST與nZST在c=0.4～0.5時(shí)功效類似，而且二者始終優(yōu)于ZST．這說明即使β值相關(guān)性較弱時(shí)，我們提出的協(xié)方差修正方法仍有優(yōu)秀的表現(xiàn)．

圖3 3種方法第一類錯(cuò)誤率與功效隨參數(shù)c的變化Fig.3 Type I error and power comparison among three methods under different values of ccZST、ZST、nZST第一類錯(cuò)誤率(a)與功效(b)與隨參數(shù)c值變化情況(α=0.05).3種方法均能將第一類錯(cuò)誤率控制在0.05左右．c≥0.4時(shí)cZST即可作為首選方法．T0=15，其余參數(shù)同實(shí)驗(yàn)1．

實(shí)驗(yàn)3： 探究ZST、cZST和nZST在測序讀數(shù)小的數(shù)據(jù)集上的功效表現(xiàn)，c取0.3令β值相關(guān)性較弱，T0分別取15、10和5，3種方法FDR結(jié)果見圖4．由實(shí)驗(yàn)2我們已知，T0=15，c=0.3時(shí)，nZST功效優(yōu)于cZST，能將FDR控制在更低的范圍，如圖4左圖所示．但當(dāng)T0減小時(shí)(圖4中、右)，同水平下cZST的FDR變得比nZST更低，意味著當(dāng)區(qū)域測序讀數(shù)降低時(shí)，cZST方法將具有格外的優(yōu)勢．

圖4 3種方法在測序總讀數(shù)不同時(shí)FDR比較Fig.4 FDR values of three methods under different reading depths模型的泊松分布期望T0變化時(shí)cZST、ZST、nZST方法的FDR值．T0=15(左)，10(中)，5(右)．c=0.3，其余參數(shù)同實(shí)驗(yàn)1．

實(shí)驗(yàn)4： 將cZST、nZST與已有方法比較．我們提出的方法能夠以基因功能區(qū)域?yàn)閱挝?，量化甲基化水平差異程度，獲得單一p值，供疾病基因關(guān)聯(lián)研究．目前，適用BS-seq數(shù)據(jù)，在染色體區(qū)域上搜索差異甲基化的方法有COHCAP[10]、BISEQ[11]、HMM-Fisher[12]等．它們使用的統(tǒng)計(jì)方法、拼接區(qū)域的標(biāo)準(zhǔn)不盡相同，方法之間難以互相比較．因此我們采用Dolzhenko等[17]在BS-seq研究中使用的Stouffer-Liptak檢驗(yàn)[18]，計(jì)算單位點(diǎn)檢驗(yàn)之間的協(xié)方差，通過合并已有方法在基因功能區(qū)域內(nèi)平行的多個(gè)DML檢驗(yàn)p值，最終得到單一p值，與我們的方法進(jìn)行功效比較．

我們將功效表現(xiàn)較好的cZST、nZST與上面3種已有方法進(jìn)行功效比較．區(qū)域位點(diǎn)總個(gè)數(shù)L從4變化到16，令μ1=μ0=0產(chǎn)生數(shù)據(jù)，計(jì)算第一類錯(cuò)誤率，再調(diào)整μ1中央的兩個(gè)值為1，即令區(qū)域中央2個(gè)位點(diǎn)甲基化水平存在差異，比較各方法功效，結(jié)果如圖5所示．在能夠控制住第一類錯(cuò)誤率的前提下，cZST檢驗(yàn)功效隨L增大而緩落，在位點(diǎn)總數(shù)超過8個(gè)以后表現(xiàn)優(yōu)于其他所有方法．說明cZST方法是一種靈敏度較高的方法，功效主要取決區(qū)域內(nèi)存在差異甲基化位點(diǎn)的個(gè)數(shù)和差異情況，受無關(guān)干擾位點(diǎn)個(gè)數(shù)影響較?。?/p>

圖5 5種方法第一類錯(cuò)誤率與統(tǒng)計(jì)功效隨區(qū)域內(nèi)位點(diǎn)個(gè)數(shù)變化的趨勢Fig.5 Type I error rate and power comparison among methods under different site numbers in a region存在甲基化差異的位點(diǎn)個(gè)數(shù)設(shè)為0(a)或2(b)，令L從4變化到16產(chǎn)生模擬數(shù)據(jù)，分別比較5種方法第一類錯(cuò)誤率與功效的變化情況(α=0.05)．c=1，T0=15，其他參數(shù)同實(shí)驗(yàn)1．

綜上，針對基因功能區(qū)域的cZST檢驗(yàn)有著較高的靈敏度，當(dāng)區(qū)域內(nèi)擾動(dòng)位點(diǎn)多時(shí)，功效優(yōu)于已有方法．且在兩位點(diǎn)β值的相關(guān)系數(shù)與距離有關(guān)的情況下，無論測序平均讀數(shù)多少，該種策略相比原始ZST，均能在控制住第一類錯(cuò)誤率的同時(shí)帶來功效提升．分析測序平均總讀數(shù)低的數(shù)據(jù)時(shí)，檢驗(yàn)優(yōu)勢更為明顯．

3 真實(shí)數(shù)據(jù)分析

我們選用來自一項(xiàng)乳腺導(dǎo)管癌研究的GEO編號為GSE69993[19]的甲基化RRBS數(shù)據(jù)，應(yīng)用cZST進(jìn)行差異甲基化檢驗(yàn)．考慮到環(huán)境因素對DNA甲基化測序的影響，我們使用同一批次的15個(gè)病例數(shù)據(jù)和5個(gè)對照數(shù)據(jù)．首先，對415985個(gè)候選位點(diǎn)進(jìn)行預(yù)處理，丟棄平均讀數(shù)少于4個(gè)的位點(diǎn)．其次，對剩余位點(diǎn)使用TxDb.Hsapiens.UCSC.hg19.knownGene程序注釋基因信息，包括基因名稱以及是否處于啟動(dòng)子、轉(zhuǎn)錄起始位點(diǎn)、外顯子及轉(zhuǎn)錄終止位點(diǎn)上．利用注釋信息我們將位點(diǎn)劃歸區(qū)域，注釋信息與前一個(gè)位點(diǎn)不同或間距超過1000bp的CpG位點(diǎn)將被劃入下一個(gè)新區(qū)域．此外，為保證檢驗(yàn)的可靠性及運(yùn)算速度，將位點(diǎn)數(shù)超過40的區(qū)域平均分為兩個(gè)小區(qū)域．對所有區(qū)域進(jìn)行標(biāo)號，包含多位點(diǎn)區(qū)域的數(shù)目與包含單位點(diǎn)區(qū)域的數(shù)目比值約為0.93．最終我們篩選出32360個(gè)區(qū)域，修正協(xié)方差矩陣，并進(jìn)行多位點(diǎn)得分檢驗(yàn)；在對篩選出的34682個(gè)單位點(diǎn)使用得分檢驗(yàn)時(shí)，由于只有一個(gè)位點(diǎn)，cZST與nZST是等價(jià)的．最后進(jìn)行104次置換檢驗(yàn)，來計(jì)算p值．同時(shí)，運(yùn)用COHCAP對這些區(qū)域和單位點(diǎn)進(jìn)行差異甲基化檢驗(yàn)，并與cZST結(jié)果進(jìn)行比較．

圖6展示了cZST方法的運(yùn)行結(jié)果，圖上標(biāo)示出了5個(gè)乳腺癌關(guān)聯(lián)基因的p值水平，這5個(gè)基因均出現(xiàn)在Cancer Genetics Web乳腺癌基因集中，并有大量研究表明其與乳腺癌的關(guān)聯(lián)性[20-24]，這說明cZST能有效篩選出疾病關(guān)聯(lián)的基因功能區(qū)域．當(dāng)p值的閾值為0.0001時(shí)，則cZST檢測出672個(gè)區(qū)域和129個(gè)單位點(diǎn)，涉及3560個(gè)CpG位點(diǎn)．隨后把間距小于1000bp且注釋相同的區(qū)域和單位點(diǎn)進(jìn)行合并，得到661個(gè)基因功能區(qū)和81個(gè)單位點(diǎn)，合計(jì)742個(gè)顯著結(jié)果．基于同樣閾值下，COHCAP檢測出2920個(gè)顯著性CpG位點(diǎn)，利用基因注釋信息分為571個(gè)區(qū)域和168個(gè)單位點(diǎn)共739個(gè)顯著結(jié)果．cZST顯著性功能區(qū)域與單位點(diǎn)個(gè)數(shù)比值為8.16，遠(yuǎn)高于待檢時(shí)的0.93，也高于COHCAP的3.40，說明cZST方法在基因功能區(qū)域方面更具優(yōu)勢．

圖6 cZST檢驗(yàn)結(jié)果的曼哈頓圖Fig.6 Manhattan plot of cZST test results檢驗(yàn)結(jié)果顯著性水平在各染色體中的分布情況．各位點(diǎn)圖中標(biāo)注了5個(gè)乳腺癌關(guān)聯(lián)基因，其基因功能區(qū)檢驗(yàn)p值均≤10-4．

從檢驗(yàn)顯著性結(jié)果相對CpG島的位置分類來看，cZST方法76.3%的顯著性位點(diǎn)位于CpG島上，COHCAP則為77.5%，兩種方法檢驗(yàn)結(jié)果都說明大部分甲基化對表達(dá)調(diào)控是通過CpG島的甲基化或去甲基化實(shí)現(xiàn)的，這一點(diǎn)符合以往研究結(jié)論[25]．顯著性結(jié)果基因注釋信息如圖7所示，cZST方法約43.3%的顯著性結(jié)果能被hg19基因注釋覆蓋(維恩圖彩色部分)，此比例在COHCAP方法中約為40.7%．在外顯子和轉(zhuǎn)錄起始位點(diǎn)占比上，cZST結(jié)果基本與COHCAP的結(jié)果持平，啟動(dòng)子占比cZST遠(yuǎn)高于COHCAP．結(jié)合以往研究結(jié)論，啟動(dòng)子或轉(zhuǎn)錄起始位點(diǎn)發(fā)生甲基化對基因表達(dá)影響更大，因此甲基化水平差異也應(yīng)大概率地出現(xiàn)在這兩個(gè)區(qū)域，這從側(cè)面說明cZST結(jié)果比COHCAP更為可靠．

圖7 cZST與COHCAP顯著性結(jié)果注釋信息維恩圖Fig.7 Annotation Venn diagrams of cZST and COHCAP results圖中空白部分?jǐn)?shù)字代表沒有基因注釋標(biāo)記的顯著性結(jié)果占比．

將cZST與COHCAP檢驗(yàn)p值從小到大排序，取前2000個(gè)基因功能區(qū)域，使用DAVID工具進(jìn)行富集分析并聚類功能注釋[26]，按相似程度整合得到若干功能大類．選取兩種方法前4個(gè)的功能注釋聚類，按顯著性水平將注釋詞條排序，結(jié)果見表1、表2．

表1 cZST顯著結(jié)果基因富集分析功能注釋集Tab.1 Functional annotation clusters of gene enrichment analysis with cZST results

注： 表中E-x表示×10-x．

表2 COHCAP顯著結(jié)果基因富集分析功能注釋集Tab.2 Functional annotation clusters of gene enrichment analysis with COHCAP results

(續(xù)表)

注： 表中E-x表示×10-x．

cZST對應(yīng)富集分析結(jié)果前4類基因注釋集依次為： PH結(jié)構(gòu)域、基因轉(zhuǎn)錄調(diào)控、Rho蛋白DH結(jié)構(gòu)域以及鈣黏著蛋白；COHCAP對應(yīng)富集分析結(jié)果前4類基因注釋集依次為： 基因轉(zhuǎn)錄調(diào)控、鈣黏著蛋白、ATP結(jié)合以及細(xì)胞黏連．方法在功能注釋集： 基因轉(zhuǎn)錄調(diào)控和鈣黏著蛋白詞條上有著類似結(jié)果，可以相互驗(yàn)證．基因轉(zhuǎn)錄調(diào)控毋庸置疑是DNA甲基化最重要的生物功能．文獻(xiàn)資料表明鈣黏著蛋白的表達(dá)直接影響腫瘤的浸潤和轉(zhuǎn)移[27]，已有研究證實(shí)了E-黏著蛋白啟動(dòng)子甲基化調(diào)節(jié)該蛋白表達(dá)水平，并與乳腺導(dǎo)管癌患病狀態(tài)及生存時(shí)間關(guān)聯(lián)[28-29]．以上兩個(gè)功能聚類的結(jié)果在cZST富集分析中排名第二、第四，這佐證了該方法的效用．

而cZST對應(yīng)的第一、第三類富集結(jié)果指向了小G蛋白Rho調(diào)控及影響小G蛋白構(gòu)型的PH結(jié)構(gòu)域(Pleckstrin homolgy domain)、DH結(jié)構(gòu)域(Dbl homolgy domain)相關(guān)基因．已有很多研究探討了在乳腺導(dǎo)管癌樣本中Rho表達(dá)與腫瘤轉(zhuǎn)移狀況的相關(guān)性和可能的生物學(xué)機(jī)理[30-31]．與COHCAP對應(yīng)結(jié)果中第三注釋集ATP結(jié)合功能、第四類注釋集細(xì)胞黏連相比，cZST的基因富集功能聚類與乳腺癌關(guān)聯(lián)更大．并且cZST顯著性基因富集分析結(jié)果說明，Rho調(diào)控基因、Pleckstrin及Dbl基因家族的基因功能區(qū)域甲基化，可能影響小G蛋白Rho的構(gòu)型和活性，進(jìn)而在乳腺導(dǎo)管癌的致病或擴(kuò)散惡化過程中起作用，對后續(xù)研究起到較好的指引作用．

4 總結(jié)與展望

本研究中，我們從觀察真實(shí)數(shù)據(jù)出發(fā)，基于染色體上同間距甲基化位點(diǎn)的相關(guān)系數(shù)離散度小這一特征，提出了利用位點(diǎn)間物理距離，引入甲基化水平相關(guān)系數(shù)的全局信息，對協(xié)方差進(jìn)行修正的方法cZST．該方法在研究低讀數(shù)、無差異位點(diǎn)多的BS-seq數(shù)據(jù)時(shí)展示出優(yōu)于已有方法的功效．通過模擬分析和真實(shí)RRBS數(shù)據(jù)研究，可以看出cZST方法在基因功能區(qū)域甲基化差異分析中具有較高靈敏度，并且適合分析測序平均總讀數(shù)較低的數(shù)據(jù)．在真實(shí)數(shù)據(jù)處理分析中，對顯著性結(jié)果進(jìn)行富集分析，cZST的基因集對應(yīng)的p值顯著功能注釋集，擁有與乳腺導(dǎo)管癌更緊密的聯(lián)系．并且作為一項(xiàng)前瞻性研究，找到了調(diào)控Rho基因及Pleckstrin、Dbl基因家族的甲基化與該疾病的密切聯(lián)系，可為后續(xù)設(shè)計(jì)實(shí)驗(yàn)探究乳腺導(dǎo)管癌遺傳機(jī)理提供幫助．

本方法仍有幾個(gè)方面值得改進(jìn)： (1) cZST修正方法局限于協(xié)方差矩陣，未能采用如BSmooth方法[32]的策略，利用臨近位點(diǎn)甲基化信息對讀數(shù)小的β值進(jìn)行修正.(2) 此方法在處理區(qū)域內(nèi)位點(diǎn)間距遠(yuǎn)、相關(guān)系數(shù)較小的稀疏數(shù)據(jù)時(shí)，比如模擬實(shí)驗(yàn)中c<0.4的情況，表現(xiàn)功效可能不如nZST．因此應(yīng)根據(jù)數(shù)據(jù)特征決定選用方法.(3) 本方法計(jì)算p值時(shí)采用置換檢驗(yàn)，在基因功能區(qū)域協(xié)方差、β值確定的情況下，檢驗(yàn)的時(shí)間取決于置換次數(shù)．因此如果多重檢驗(yàn)需要對p值作Bonferroni校正時(shí)必須保證進(jìn)行相當(dāng)數(shù)量級的置換次數(shù)，檢驗(yàn)速度將受到一定影響．

5 附 錄

將全部樣本甲基化數(shù)據(jù)β分為病例組βA(affected)與對照組βN(normal)．得分向量U的等價(jià)表達(dá)式為：

從而有

[1] MCCABE M T, DAVIS J N, DAY M L. Regulation of DNA methyltransferase1 by the pRb/E2F1 pathway [J].CancerResearch, 2005,65(9)： 3624-3632.

[2] VARRIALE A, BERNARDI G. Distribution of DNA methylation, CpGs, and CpG islands in human isochores [J].Genomics, 2010,95(1)： 25-28.

[3] MA X, WANG Y W, ZHANG M Q,etal. DNA methylation data analysis and its application to cancer research [J].Epigenomics, 2013,5(3)： 301-316.

[4] KISSELJOVA N P, ZUEVA E S, PEVZNER V S,etal. De novo methylation of selective CpG dinucleotide clusters in transformed cells mediated by an activated N-ras [J].InternationalJournalofOncology, 1998,12： 203-209.

[5] JONES P A. Functions of DNA methylation： islands, start sites, gene bodies and beyond [J].NatureReviewsGenetics, 2012,13(7)： 484-492.

[6] SKINNER M K, GUERRERO-BOSAGNA C. Role of CpG deserts in the epigenetic transgenerational inheritance of differential DNA methylation regions [J].BMCGenomics, 2014,15(1)： 692.

[7] LI N, YE M, LI Y,etal. Whole genome DNA methylation analysis based on high throughput sequencing technology [J].Methods, 2010,52(3)： 203-212.

[8] MEISSNER A, GNIRKE A, BELL G W,etal. Reduced representation bisulfite sequencing for comparative high-resolution DNA methylation analysis [J].NucleicAcidsResearch, 2005,33(18)： 5868-5877.

[9] ANIRUDDHA C, RODGER E J, STOCKWELL P A,etal. Technical considerations for reduced representation bisulfite sequencing with multiplexed libraries [J].JournalofBiomedicine&Biotechnology, 2012,2012(4)： 741542.

[10] WARDEN C D, LEE H, TOMPKINS J D,etal. COHCAP： An integrative genomic pipeline for single-nucleotide resolution DNA methylation analysis [J].NucleicAcidsResearch, 2013,41(11)： e117.

[11] HEBESTREIT K, DUGAS M, KLEIN H U. Detection of significantly differentially methylated regions in targeted bisulfite sequencing data [J].Bioinformatics, 2013,29(13)： 1647-1653.

[12] SUN S, YU X. HMM-Fisher： identifying differential methylation using a hidden Markov model and Fisher’s exact test [J].StatisticalApplicationsinGenetics&MolecularBiology, 2016,15(1)： 55-67.

[13] ROBINSON M D, KAHRAMAN A, LAW C W,etal. Statistical methods for detecting differentially methylated loci and regions [J].FrontiersinGenetics, 2014,5： 324.

[14] PAN W, KIM J, ZHANG Y,etal. A powerful and adaptive association test for rare variants [J].Genetics, 2014,197(4)： 1081-1095.

[15] SCHOOFS T, ROHDE C, HEBESTREIT K,etal. DNA methylation changes are a late event in acute promyelocytic leukemia and coincide with loss of transcription factor binding [J].Blood, 2013,121(1)： 178-187.

[16] WANG C, SUN L, ZHENG H,etal. Detecting multiple variants associated with disease based on sequencing data of case-parent trios [J].JournalofHumanGenetics, 2016,61(10)： 851-860.

[17] DOLZHENKO E, SMITH A D. Using beta-binomial regression for high-precision differential methylation analysis in multifactor whole-genome bisulfite sequencing experiments [J].BMCBioinformatics, 2014,15(1)： 215.

[19] ABBA M C, GONG T, LU Y,etal. A molecular portrait of high-grade ductal carcinomainsitu[J].CancerResearch, 2015,75(18)： 3980-3990.

[20] SLATTERY M L, LUNDGREEN A, JOHN E M,etal. MAPK genes interact with diet and lifestyle factors to alter risk of breast cancer： The breast cancer health disparities study [J].Nutrition&Cancer, 2015,67(2)： 1-13.

[21] GRAHAM T R, YACOUB R, TALIAFERROSMITH L,etal. Reciprocal regulation of ZEB1 and AR in triple negative breast cancer cells [J].BreastCancerResearchandTreatment, 2010,123(1)： 139-147.

[22] SILBERSTEIN G B, VAN H K, STRICKLAND P,etal. Altered expression of the WT1 wilms tumor suppressor gene in human breast cancer [J].ProcNatlAcadSciUSA, 1997,94(15)： 8132-8137.

[23] CHENG Q, CHANG J T, GERADTS J,etal. Amplification and high-level expression of heat shock protein 90 marks aggressive phenotypes of human epidermal growth factor receptor 2 negative breast cancer [J].BreastCancerResearch, 2012,14(2)： R62.

[24] PEDERSON H J, PADIA S A, MAY M,etal. Managing patients at genetic risk of breast cancer [J].CleveClinJMed, 2016,83(3)： 199-206.

[25] YONG W, LEUNG F C C. An evaluation of new criteria for CpG islands in the human genome as gene markers [J].Bioinformatics, 2004,20(7)： 1170-1177.

[26] HUANG D W, SHERMAN B T, LEMPICKI R A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources [J].NatureProtocol, 2009,4(1)： 44-57.

[27] FRIXEN U H, BEHRENS J, SACHS M,etal. E-cadherin-mediated cell-cell adhesion prevents invasiveness of human carcinoma cells [J].JournalofCellBiology, 1991,113(1)： 173-185.

[28] HU X C, LOO W T, CHOW L W. E-cadherin promoter methylation can regulate its expression in invasive ductal breast cancer tissue in Chinese woman [J].LifeSciences, 2002,71(12)： 1397-1404.

[30] HALON A, DONIZY P, SUROWIAK P,etal. ERM/Rho protein expression in ductal breast cancer： A15 year follow-up [J].CellularOncology, 2013,36(3)： 181-190.

[31] MURAKAMI E, NAKANISHI Y, HIROTANI Y,etal. Roles of ras homolog a in invasive ductal breast carcinoma [J].ActaHistochemicaEtCytochemica, 2016,49(5)： 131-140.

[32] HANSEN K D, IRIZARRY R A, ZHIJIN W. Removing technical variability in RNA-seq data using conditional quantile normalization [J].Biostatistics, 2012,13(2)： 204-216.

AStatisticalMethodforDetectingDifferentiallyMethylatedGeneFunctionalRegionsBasedonCase-ControlStudies

MAXinyu,WANGChan,HUYueqing

(DepartmentofBiostatisticsandComputationalBiology,SchoolofLifeSciences,FudanUniversity,Shanghai200433,China)

DNA methylation is an important epigenetic modification that participates in the regulation of gene expression and other biological activities. In case-control studies, searching for differentially methylated gene functional regions associated with disease can provide clues of disease-related gene list for follow-up studies. Nowadays existing differential methylation detection methods are lack of a uniform standard for defining differentially methylated regions. Also test results can be across multiple genes or multiple functional regions, bringing difficulty for biological interpretation. In this paper, a statistical approach, cZST, is proposed based on Z-score test. It is able to detect differences of methylation levels in each gene functional region. This method uses physical distances between sites and the global information of methylation levels in each chromosome to modify the covariance matrix. The simulation results show that cZST is superior to other methods in BS-seq data of a small sample size or a low reading depth. In the real data analysis of a breast ductal cancer study, 661 significant functional regions and 81 significant single loci are recognized. Gene annotation analysis of results demonstrates the effectiveness of cZST. Also we find that in regulator genes of Rho,Dbl and Pleckstrin gene families, DNA methylation level differences are associated with breast ductal carcinoma affection.

differential methylation detection; BS-seq; functional enrichment analysis

0427-7104(2017)06-0701-11

2017-04-05

國家自然科學(xué)基金(11571082，11171075)

馬欣宇(1992—)，男，碩士研究生；胡躍清，男，教授，通信聯(lián)系人，E-mail： yuehu@fudan.edu.cn

Q348

A