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

馬欣宇，王 嬋，胡躍清

(復旦大學 生命科學學院 生物統(tǒng)計學與計算生物學系，上海 200433)

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

馬欣宇，王 嬋，胡躍清

(復旦大學 生命科學學院 生物統(tǒng)計學與計算生物學系，上海 200433)

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

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

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

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

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

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

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

1 符號與方法

(1)

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

(2)

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

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

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

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

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

(3)

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

(4)

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

2 隨機模擬研究

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

σil=exp(-dil/d0),

(5)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(續(xù)表)

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

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

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

4 總結與展望

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

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

5 附 錄

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

從而有

[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

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

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

Q348

A