表观基因组甲基化数据的统计分析方法*

南京医科大学生物统计学系(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)；江苏省青蓝工程资助项目