表觀基因組甲基化數據的統計分析方法*

南京醫科大學生物統計學系(211166)

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

?

·綜述·

表觀基因組甲基化數據的統計分析方法*

南京醫科大學生物統計學系(211166)

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

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

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

本文對近幾年國內外學者提出的甲基化數據統計分析新方法進行了綜述和探討，將這些方法細分為單位點的關聯研究和多位點的關聯研究兩大類。系統介紹了每一種方法的背景、基本思想和優缺點。

單位點的分析方法

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

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

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

2.考慮年齡協變量的方法

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

(1)參數法

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

(2)非參數法

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

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

3.廣義指數傾斜模型半參數檢驗

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

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

4.Bayesian分層模型

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

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

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

多位點的分析方法

1.懲罰logistic回歸

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

懲罰logistic回歸是對Li等[26-27]提出的graph-constrained過程的改進。對位點間相關性的懲罰有兩種形式：環狀網絡和全關聯網絡。當基因內部的CpG位點之間存在相關性時，懲罰logistic回歸要優于現有的主流正則化模型，如lasso[28]、Enet[29]。位點間相關結構的選取以及該法的優劣取決于基因內CpG位點的潛在真實相關性，而這種相關性是不固定的，并且要比上文所假設的兩種結構復雜得多。此外，一個通路上的基因之間可能也存在相關性，懲罰logistic回歸未能考慮這一問題。

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

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

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

3.空間聚類法

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

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

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

存在的問題及展望

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

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

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

[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.

(責任編輯：郭海強)

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

*:國家自然基金(No.81530088,81473070,81373102,61301251,81402764)；公共衛生與預防醫學江蘇省高校優勢學科建設專項資金資助；江蘇省高等學校自然科學項目(No.12KJB310003)；江蘇省青藍工程資助項目