表觀基因組甲基化數(shù)據(jù)的統(tǒng)計分析方法*

南京醫(yī)科大學(xué)生物統(tǒng)計學(xué)系(211166)

張秋伊　陳　峰　魏永越　郭　麗　趙　楊△

?

·綜述·

表觀基因組甲基化數(shù)據(jù)的統(tǒng)計分析方法*

南京醫(yī)科大學(xué)生物統(tǒng)計學(xué)系(211166)

張秋伊陳峰魏永越郭麗趙楊△

近年來，全基因組關(guān)聯(lián)研究(genome-wide association studies,GWAS)取得了舉世矚目的成就，識別了成千上萬個復(fù)雜疾病(complex disease)遺傳易感位點[1]。然而，疾病的發(fā)生還與遺傳以外的諸多因素有關(guān)，包括表觀遺傳(epigenetics)的改變[2]。這類問題的闡述需要通過與全基因組關(guān)聯(lián)研究類似的大規(guī)模、系統(tǒng)性的研究——表觀基因組關(guān)聯(lián)研究(epigenome-wide association studies,EWAS)[3-4]。其中，DNA甲基化(DNA methylation)與人類發(fā)育和腫瘤疾病密切相關(guān)，已成為表觀遺傳學(xué)的主要研究內(nèi)容[5-6]。哺乳動物的DNA甲基化通常發(fā)生在DNA鏈中胞嘧啶-鳥嘌呤核苷酸對(CpGs)的胞嘧啶上，生成5-甲基胞嘧啶。其改變的是基因表達(dá)水平，而不是DNA序列，這種改變是可遺傳的。它能記錄人體一生的環(huán)境暴露，為疾病診斷和危險因素分層提供有用的生物標(biāo)志信息來源。

新一代測序和生物芯片技術(shù)的發(fā)展為表觀基因組研究創(chuàng)造機(jī)遇的同時，又給數(shù)據(jù)處理、統(tǒng)計分析和結(jié)果的生物學(xué)解釋帶來很大挑戰(zhàn)[7-8]。甲基化數(shù)據(jù)分析的主要目的是識別差異甲基化位點(differentially methylated loci,DML)。甲基化水平代表值主要有β-值和M值[9]，β-值為0～1的比例，從生物學(xué)角度易于理解，但通常每個值間方差不齊；β-值經(jīng)過對數(shù)變換為M值，可達(dá)到方差齊性，但方差不齊可能蘊含的生物學(xué)信息，如批次效應(yīng)(batch effect)，會在變換過程中損失[10-11]。目前常用的統(tǒng)計分析方法主要有t檢驗、基于回歸的方法和基于秩次的方法[12]。EWAS研究亟需DNA甲基化數(shù)據(jù)的標(biāo)準(zhǔn)化算法和識別病例對照差異甲基化位點的穩(wěn)健統(tǒng)計分析方法[13]。

本文對近幾年國內(nèi)外學(xué)者提出的甲基化數(shù)據(jù)統(tǒng)計分析新方法進(jìn)行了綜述和探討，將這些方法細(xì)分為單位點的關(guān)聯(lián)研究和多位點的關(guān)聯(lián)研究兩大類。系統(tǒng)介紹了每一種方法的背景、基本思想和優(yōu)缺點。

單位點的分析方法

1.基于均勻-正態(tài)混合分布模型的似然比檢驗

來自Illumina芯片的β-值通常呈現(xiàn)雙峰分布，峰值位于完全甲基化(βj=1)和未甲基化(βj=0)處。盡管DNA甲基化水平值是0～1的定量資料，但在分子水平，位點的甲基化狀態(tài)有未甲基化、完全甲基化和半甲基化(只有一側(cè)的胞嘧啶甲基化)三種[14]。因此，Wang[15]于2011年提出了基于均勻-正態(tài)混合分布模型的似然比檢驗方法，來識別病例和對照間的差異甲基化位點。其基本思想：根據(jù)甲基化位點的三種不同狀態(tài)，通過三組分混合分布(兩均勻分布和一截斷的正態(tài)分布)來模擬甲基化數(shù)據(jù)。通過混合分布的概率和正態(tài)分布的均數(shù)來檢驗差異甲基化位點。

當(dāng)病例組和對照組整體的甲基化水平均數(shù)接近，混合概率和正態(tài)分布均數(shù)存在差別時，該方法優(yōu)于t檢驗。但這種方法過于保守，并且EM算法迭代過程運行速度慢，需要占用較多計算資源。

2.考慮年齡協(xié)變量的方法

研究表明甲基化水平與年齡存在著很大關(guān)系[16-17]，為調(diào)整年齡這一混雜因素，Chen等提出了幾種解決方法[18-20]。

(1)參數(shù)法

基本思想：首先將樣本分為若干個年齡組，每個年齡分組對病例組和對照組進(jìn)行兩次方差不等的t檢驗，獲得兩次單側(cè)檢驗的P值，根據(jù)這兩組P值估計該位點總的P值[18]。假設(shè)共有k個年齡組，左側(cè)檢驗的P值用Pli(i=1,2,…,k)表示，相應(yīng)的右側(cè)檢驗的P值用Pri表示。根據(jù)Fisher合并檢驗[21]，可得到：

(2)非參數(shù)法

考慮到β-值分布的非正態(tài)性，作者提出非參數(shù)法來代替上文的t檢驗[19]。基本思想：同樣將樣本分為若干個年齡組，每個年齡分組對病例組和對照組進(jìn)行非參數(shù)Kruskal-Wallis(KW)檢驗，獲得每個年齡組比較的P值后估計得到整體的P值。整體P值的估計仍采用Fisher合并檢驗，該方法也可用于多個樣本資料的比較(如對照組、療前組、療后組)。

考慮到多個樣本間甲基化水平存在增大或減小的趨勢，例如在對照、療后和療前這三組間甲基化水平逐漸降低，Chen等[20]提出將Cuzick非參數(shù)趨勢性檢驗用于此類數(shù)據(jù)，得到單側(cè)檢驗P值后估計整體P值，此方法可獲得較高的檢驗效能。

3.廣義指數(shù)傾斜模型半?yún)?shù)檢驗

有研究發(fā)現(xiàn)，不同組間甲基化水平的方差也存在差異[8,21]，方差不齊可能蘊含批次效應(yīng)等生物學(xué)信息，在統(tǒng)計分析過程中需要保留這些信息，因此Chen等[22]于2013年提出精簡的兩樣本廣義指數(shù)傾斜模型。該方法為半?yún)?shù)方法，首先假設(shè)兩組甲基化數(shù)據(jù)服從相同分布，建立比較模型，來捕獲均數(shù)和方差之間的差別[23]。

相對于t檢驗和基于回歸的方法而言，該法還可以識別兩組數(shù)據(jù)方差的差別；而當(dāng)患者與正常人甲基化水平僅存在均數(shù)差別時，該法的檢驗效能低于t檢驗。基于指數(shù)傾斜模型的經(jīng)驗對數(shù)似然比檢驗和偽似然比檢驗可以利用方差不齊所包含的信息，作為傳統(tǒng)方法的補充。

4.Bayesian分層模型

Feng等[24](2014)提出Bayesian分層模型的方法，采用beta-二項分布分層模型來解決不同分組CpG位點甲基化水平方差不齊和樣本量較小的問題。

nφ(φij-1-1)(1-μij)+nφ(φij-1-1)-

分層模型中，beta分布用于解釋個體間的生物學(xué)變異，二項分布則解釋測序過程中DNA片段隨機(jī)抽樣帶來的測量誤差。尤其小樣本情況下，這種方法明顯優(yōu)于其他常用方法。除病例對照研究外，Beyesian分層模型也可用于更復(fù)雜的試驗設(shè)計，如多組比較、連續(xù)性結(jié)局變量等。

多位點的分析方法

1.懲罰logistic回歸

一個基因內(nèi)的CpG位點的甲基化水平通常存在相關(guān)性。這些相關(guān)的位點中，一部分是致病位點，而另外一些位點是中性的。根據(jù)這一特點，Sun等[25](2012)提出懲罰logistic回歸模型來篩選基因內(nèi)相關(guān)的CpG位點。這種方法在考慮相關(guān)性的前提下，獨立篩選疾病相關(guān)CpG位點。

懲罰logistic回歸是對Li等[26-27]提出的graph-constrained過程的改進(jìn)。對位點間相關(guān)性的懲罰有兩種形式：環(huán)狀網(wǎng)絡(luò)和全關(guān)聯(lián)網(wǎng)絡(luò)。當(dāng)基因內(nèi)部的CpG位點之間存在相關(guān)性時，懲罰logistic回歸要優(yōu)于現(xiàn)有的主流正則化模型，如lasso[28]、Enet[29]。位點間相關(guān)結(jié)構(gòu)的選取以及該法的優(yōu)劣取決于基因內(nèi)CpG位點的潛在真實相關(guān)性，而這種相關(guān)性是不固定的，并且要比上文所假設(shè)的兩種結(jié)構(gòu)復(fù)雜得多。此外，一個通路上的基因之間可能也存在相關(guān)性，懲罰logistic回歸未能考慮這一問題。

2.高分辨率甲基化譜的整體分析

除了差異甲基化位點外，有時我們也關(guān)注整個表觀基因組的甲基化水平差異。例如，對于癌癥和年齡相關(guān)疾病呈現(xiàn)的是整個基因組DNA的低甲基化狀態(tài)。因此，Zhao等[30]于2015年提出了針對表觀基因組或者許多個CpG位點甲基化譜的整體分析方法(global analysis of methylation profiles,GAMP)。

其原理是整體甲基化差別體現(xiàn)在CpG甲基化水平整體分布的差異，少數(shù)位點甲基化水平的改變不會對整個分布產(chǎn)生很大影響。用B-Spline系數(shù)來概括甲基化值的整體分布，采用方差成分檢驗整體甲基化水平的差別。兩組間系數(shù)差別的檢驗采用方差成分檢驗[31-32]。其優(yōu)點在于自由度取決于回歸模型系數(shù)間的相關(guān)性，若相關(guān)性高，則自由度較小，從而提高檢驗效能；另一方面還可以將需要調(diào)整的協(xié)變量納入回歸模型。該方法可用于整個表觀基因組甲基化的整體分析，此外，為方便結(jié)果的解釋，也可將CpG位點限制于相關(guān)功能區(qū)域，包括如CpG島、啟動子區(qū)等。但這種方法適用于檢驗整體甲基化水平的差異，若位點數(shù)很少，就不足以估計概率密度和CDF，因此，作者要求CpG位點數(shù)達(dá)到50以上。

3.空間聚類法

有學(xué)者指出甲基化水平是叢集的，如啟動子區(qū)的甲基化位點共同影響基因表達(dá)水平[33]。利用位點間距離的信息，在關(guān)聯(lián)研究中我們就可以獲得更高的檢驗效能。Yip等[34](2014)提出空間聚類法(spatial clustering method,SCM)，來尋找基因組中與疾病有關(guān)的候選差異甲基化區(qū)域。

空間數(shù)據(jù)分析要求資料包含區(qū)域信息，即每個位點的位置和位點間的距離。CpG位點可看成當(dāng)染色體被拉直后，沿著染色體排列的點。通過芯片測序數(shù)據(jù)，可以得到單個位點的甲基化值，該方法需要將這些甲基化值轉(zhuǎn)化為甲基化單位。對每個位點的轉(zhuǎn)換需用到一個權(quán)重：位點間距離越接近，甲基化水平越低，權(quán)重就越高。該權(quán)重既考慮了位點間距離越近，甲基化水平相關(guān)性越高的特點，又調(diào)整了位點間甲基化水平的不均勻性。分別對病例和對照組計算距離向量，表示甲基化單位的距離分布。零假設(shè)為兩組的距離分布相同，采用Ansari-Bradley非參數(shù)檢驗。

SCM在構(gòu)建統(tǒng)計量時，既包括了位點的甲基化值，又包括了空間位置信息。設(shè)定包含固定CpG位點數(shù)的基因窗，從染色體起始處滑動至末尾，篩選有意義的區(qū)域，便于進(jìn)一步的分析。但檢驗統(tǒng)計量的分布要采用permutation獲得，需要消耗更多的計算資源。協(xié)變量的調(diào)整不如GAMP法方便，只能通過分層分析、匹配或者傾向性得分的方法。此外，SCM還要求數(shù)據(jù)包含位點的位置信息，密集的Illumina Infinium 450K芯片數(shù)據(jù)提供的信息比稀疏的Illumina Infinium 27K芯片數(shù)據(jù)更為豐富。

存在的問題及展望

單位點的分析方法主要著眼于DNA甲基化水平β-值是0～1之間的定量資料，不服從正態(tài)分布且方差不齊的特點，盡可能地整合數(shù)據(jù)信息，從而提高方法的檢驗效能。由于年齡與甲基化水平間存在著高度相關(guān)，在關(guān)聯(lián)研究中，如何調(diào)整年齡這一混雜因素的影響也是這些分析方法需要考慮的問題。一些研究指出CpG位點間的甲基化水平存在著相關(guān)性，并且在不同的組織和細(xì)胞類型中均有這種相關(guān)結(jié)構(gòu)。單位點的關(guān)聯(lián)研究將每個CpG位點作為單獨的因素來分析，沒有考慮位點間的相關(guān)結(jié)構(gòu)，信息利用不充分；另外，自變量的個數(shù)遠(yuǎn)遠(yuǎn)大于樣本個數(shù)，嚴(yán)格的檢驗水準(zhǔn)校正也會帶來統(tǒng)計學(xué)效能的損失。

多位點的關(guān)聯(lián)研究彌補了單位點關(guān)聯(lián)研究的不足，利用位點間的相關(guān)結(jié)構(gòu)所提供的信息，對多個CpG位點進(jìn)行綜合來識別差異甲基化區(qū)域。雖然不能完全避免多重比較的校正，但可以大大減少多重比較的次數(shù)。在這一區(qū)域內(nèi)，既包含致病位點，也包含中性位點，將一個基因、通路、啟動子區(qū)等作為一個整體來考慮，更加符合復(fù)雜疾病的致病機(jī)制。但是，CpG位點間的相關(guān)結(jié)構(gòu)較為復(fù)雜且不固定，變量間還可能存在一階或多階的交互作用，多位點關(guān)聯(lián)研究也同樣存在不能捕獲真正致病位點的風(fēng)險。除此之外，基因型和表觀基因型間的相互關(guān)系需要我們進(jìn)行綜合分析，如何把基因多態(tài)性、DNA甲基化、基因表達(dá)等信息整合起來，這將是GWAS和EWAS統(tǒng)計分析需要進(jìn)一步探討的問題。

本文所綜述的這些統(tǒng)計分析方法都有各自的適用條件，但在相同條件下，哪種方法具有更高的檢驗效能，還需要進(jìn)一步探討。

[1]Manolio TA,Collins FS.The HapMap and Genome-Wide Association Studies in Diagnosis and Therapy.Annu Rev Med,2009,60:443-456.

[2]Petronis A.Epigenetics as a unifying principle in the aetiology of complex traits and diseases.Nature,2010,465(7299):721-727.

[3]Rakyan VK,Down TA,Balding DJ，et al.Epigenome-wide association studies for common human diseases.Nat Rev Genet,2011,12(8):529-541.

[4]Egger G,Liang GN,Aparicio A,et al.Epigenetics in human disease and prospects for epigenetic therapy.Nature,2004,429(6990):457-463.

[5]Kulis M,Esteller M.DNA Methylation and Cancer.Adv Genet,2010,70:27-56.

[6]Kulis M,Queiros AC,Beekman R,et al.Intragenic DNA methylation in transcriptional regulation,normal differentiation and cancer.Bba-Gene Regul Mech,2013,1829(11):1161-1174.

[7]Laird PW.Principles and challenges of genome-wide DNA methylation analysis.Nat Rev Genet,2010,11(3):191-203.

[8]Hansen KD,Timp W,Bravo HC，et al.Increased methylation variation in epigenetic domains across cancer types.Nat Genet,2011,43(8):768-777.

[9]Saadati M,Benner A.Statistical challenges of high-dimensional methylation data.Stat Med,2014,33(30):5347-5357.

[10]Du P,Zhang XA,Huang CC,et al.Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis.BMC Bioinformatics,2010,11：587.

[11]Leek JT,Scharpf RB,Bravo HC,et al.Tackling the widespread and critical impact of batch effects in high-throughput data.Nat Rev Genet,2010,11(10):733-739.

[12]Wang D,Yan L,Hu Q,et al.IMA:an R package for high-throughput analysis of Illumina′s 450K Infinium methylation data.Bioinformatics,2012,28(5):729-730.

[13]Bock C.Analysing and interpreting DNA methylation data.Nat Rev Genet,2012,13(10):705-719.

[14]Strachan TRA.Human Molecular Genetics.3rd.New York:Garland Science,2004.

[15]Wang S.Method to Detect Differentially Methylated Loci With Case-Control Designs Using Illumina Arrays.Genet Epidemiol,2011,35(7):686-694.

[16]Christensen BC,Houseman EA,Marsit CJ,et al.Aging and Environmental Exposures Alter Tissue-Specific DNA Methylation Dependent upon CpG Island Context.Plos Genet,2009,5(8):e1000602.

[17]Teschendorff AE,Menon U,Gentry-Maharaj A,et al.Age-dependent DNA methylation of genes that are suppressed in stem cells is a hallmark of cancer.Genome Res,2010,20(4):440-446.

[18]Chen ZX,Liu QZ,Nadarajah S.A new statistical approach to detecting differentially methylated loci for case control Illumina array methylation data.Bioinformatics,2012,28(8):1109-1113.

[19]Chen ZX,Huang HW,Liu JZ,et al.Detecting differentially methylated loci for Illumina Array methylation data based on human ovarian cancer data.BMC Med Genomics,2013,6:S9.

[20]Chen ZX,Huang HW,Liu QZ.Detecting differentially methylated loci for multiple treatments based on high-throughput methylation data.BMC Bioinformatics,2014,15:142.

[21]Fisher RA.Statistical methods for research workers.4th.Edinburgh etc.:Oliver and Boyd,1932.

[22]Gervin K,Hammero M,Akselsen HE,et al.Extensive variation and low heritability of DNA methylation identified in a twin study.Genome Res,2011,21(11):1813-1821.

[23]Chen Y,Ning Y,Hong C,et al.Semiparametric Tests for Identifying Differentially Methylated Loci With Case-Control Designs Using Illumina Arrays.Genet Epidemiol,2014,38(1):42-50.

[24]Qin J.Inferences for case-control and semiparametric two-sample density ratio models.Biometrika,1998,85(3):619-630.

[25]Feng H,Conneely KN,Wu H.A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data.Nucleic Acids Res,2014,42(8):e69.

[26]Sun H,Wang S.Penalized logistic regression for high-dimensional DNA methylation data with case-control studies.Bioinformatics,2012,28(10):1368-1375.

[27]Li CY,Li HZ.Network-constrained regularization and variable selection for analysis of genomic data.Bioinformatics,2008,24(9):1175-1182.

[28]Li CY,Li HZ.Variable selection and regression analysis for graph-structured covariates with an application to genomics.The annals of applied statistics,2010,4(3):1498-1516.

[29]Tibshirani R.Regression shrinkage and selection via the lasso:a retrospective.J Roy Stat Soc B,2011,73:273-282.

[30]Zou H,Hastie T.Regularization and variable selection via the elastic net.J Roy Stat Soc B,2005,67:301-320.

[31]Zhao N,Bell DA,Maity A,et al.Global Analysis of Methylation Profiles From High Resolution CpG Data.Genet Epidemiol,2015,39(2):53-64.

[32]Wu MC,Kraft P,Epstein MP,et al.Powerful SNP-Set Analysis for Case-Control Genome-wide Association Studies.Am J Hum Genet,2010,86(6):929-942.

[33]Wu MC,Lee S,Cai TX,et al.Rare-Variant Association Testing for Sequencing Data with the Sequence Kernel Association Test.Am J Hum Genet,2011,89(1):82-93.

[34]Hackenberg M,Barturen G,Carpena P,et al.Prediction of CpG-island function:CpG clustering vs.sliding-window methods.BMC Genomics,2010,11:327.

[35]Yip WK,Fier H,DeMeo DL,et al.A Novel Method for Detecting Association Between DNA Methylation and Diseases Using Spatial Information.Genet Epidemiol,2014,38(8):714-721.

(責(zé)任編輯：郭海強)

趙楊，E-mail:zhaoyang@njmu.edu.cn

*:國家自然基金(No.81530088,81473070,81373102,61301251,81402764)；公共衛(wèi)生與預(yù)防醫(yī)學(xué)江蘇省高校優(yōu)勢學(xué)科建設(shè)專項資金資助；江蘇省高等學(xué)校自然科學(xué)項目(No.12KJB310003)；江蘇省青藍(lán)工程資助項目