高通量DNA甲基化数据的处理和分析方法

王心宇，许颖出，刘洪波，王 芳，张 岩，苏建忠

(哈尔滨医科大学生物信息科学与技术学院，黑龙江哈尔滨150080)

DNA甲基化是重要的表观遗传学修饰之一，以往的研究表明，DNA甲基化在细胞发育和分化、调控基因表达、X染色体失活、基因沉默、疾病的发生等方面扮演着重要的角色［1-3］。在真核生物中，通常是CpG二核苷酸中胞嘧啶的第五个碳原子上发生了甲基化(5mC)，胞嘧啶甲基化也可能会发生在CHG和CHH(H是除G外的任意一种核苷酸)上。全基因组的甲基化水平呈现双峰分布，而且低甲基化的区域多数是在CpG二核苷酸聚集区域(CpG岛)［4］。以往的研究发现，位于启动子区域的高甲基化的CpG岛与基因的沉默有关，可能是因为DNA甲基化阻碍转录因子结合而直接抑制了基因转录。在过去的几十年里，由于实验技术和费用的限制，DNA甲基化的数据往往只检测了基因组的局部区域，而且是低通量的数据。

二代测序技术的发展极大地推动了表观遗传调控机制的研究，基于二代测序技术发展起来的DNA甲基化的检测技术为DNA甲基化的研究提供了大量的高通量、全基因组的DNA甲基化数据。这些高通量数据的产生使得DNA甲基化研究的重点由目标基因DNA甲基化的检测转移到了全基因组DNA甲基化高通量数据的检测、存储、处理和分析上。近几年，研究者构建了多个DNA甲基化数据库，开发了大量的DNA甲基化高通量数据的处理和分析工具，使得深入的表观遗传调控机制的研究成为可能。

1 基于二代测序技术的DNA甲基化检测技术

1.1 DNA甲基化预处理方法

甲基化后的胞嘧啶(5 mC)与普通的胞嘧啶(C)在DNA序列上并无差异，如果直接使用DNA测序，将无法区分测得的胞嘧啶C是C还是5 mC。所以检测DNA甲基化需要首先对待检测的DNA序列中胞嘧啶进行预处理，将非甲基化的胞嘧啶C与甲基化的胞嘧啶5 mC区分开来，目前的DNA甲基化预处理方式主要分为三种:

(1)限制性内切酶法(Endonuclease digestion)

限制性内切酶法是指利用甲基化限制性内切酶(HpaII，MspI和HhaI等)在各自的识别位点对甲基化的胞嘧啶有不同的敏感性来检测CpG的甲基化［5］。限制性内切酶法结合二代测序的技术有MRE-seq，MCA-seq，MSCC 和 HELP-seq。尽管限制性内切酶测序法成本低、高效，然而由于检测的CpG位点局限于酶切位点附近，基因组覆盖率低，另外还存在CpG偏好性、酶切不完全导致的假阳性等问题，使用这种方法检测DNA甲基化的研究越来越少。

(2)亲和纯化法(Affinity enrichment)

亲和纯化是利用甲基化CpG结合蛋白(MBD)或者对5mC特异的抗体来亲和提纯甲基化区域。MeDIP-seq和 MBD-seq是最常用的两种结合亲和纯化和二代测序技术的DNA甲基化检测方法。基于测序的亲和纯化法能够快速、低成本地检测全基因组范围内的甲基化水平，然而它只能获得区域的甲基化水平，特别是MeDIP-seq偏向于CpG富集的区域，分散的低密度的甲基化位点可能被识别成非甲基化区域，目前还没有能够去除掉这种偏性的生物信息学方法。

(3)重亚硫酸盐转换法(Bisulphite conversion)

重亚硫酸盐转换结合二代测序技术是目前最精准的DNA甲基化检测方法，能够检测单碱基水平的甲基化状态，被称为DNA甲基化检测的“金标准”。对基因组中未发生甲基化的胞嘧啶进行重亚硫酸盐处理，将其转换成U，经PCR扩增后变成T，重亚硫酸盐转换对甲基化的胞嘧啶不起作用。通过结合二代测序，即可绘制出单碱基分辨率的全基因组DNA甲基化图谱。目前常用的重亚硫酸盐转换结合二代测序技术的DNA甲基化检测技术有BS-seq和RRBS等。

1.2 二代测序技术

在使用DNA甲基化预处理区分出未甲基化的胞嘧啶和甲基化的胞嘧啶后，再使用二代测序技术检测DNA序列，来获取胞嘧啶上的甲基化状态。

目前二代测序技术主要分为三个平台:Roche、Illumina、SOLiD。其中每种测序平台又拥有多种系统，比如Illumina就有HiSeq、GAIIx等系统。不同的测序技术在测得的read长度、精确性、通量都有差异，适用于不同的研究目的需要。

1.3 高通量DNA甲基化检测技术进展

结合二代测序技术和DNA甲基化预处理的DNA甲基化检测方法，在近几年获得了大量的全基因组的DNA甲基化测序数据。

国外很多实验室产生了大量、精准的高通量DNA甲基化数据，例如，Lister等人于2008年检测的拟南芥全基因组甲基化谱和2009年测得的人类全基因组甲基化谱［6-7］，Stadler等人于2011年测定了小鼠胚胎干细胞和神经前体细胞的全基因组甲基化谱等［8］。国内近年来也产生了大量的高通量DNA甲基化数据，例如，2010年，中科院昆明研究所，华大基因和上海交通大学癌症表观遗传中心等九家科研机构联合测定了桑蚕的单碱基水平的DNA甲基化谱，王俊教授课题组测定的人类完全分化的血细胞的全基因组DNA甲基化谱等。这些全基因组水平的DNA甲基化数据为表观遗传调控机制的研究提供了数据资源。

2 DNA甲基化数据储存和可视化

目前研究者构建了各种各样的数据库来存储世界范围的各大实验室和科研机构产生的高通量DNA甲基化数据，便于数据的查询、下载、可视化分析及全球化的资源共享。从第一个DNA甲基化的公共数据库MethDB由Grunau等人于2001年构建以来，已有多个和DNA甲基化相关的数据库被开发，例如，NCBI的存储表观遗传修饰数据的Epigenomics，主要包括DNA甲基化、组蛋白修饰和非编码RNA等数据。PubMeth是结合文本的基因注释信息的DNA甲基化数据库。DiseaseMeth储存72种人类疾病相关的DNA甲基化的数据库，并实现了统计学分析及可视化［9］。

3 高通量DNA甲基化数据的处理和分析

结合二代测序技术和DNA甲基化预处理的方法，在近几年产生了大量的全基因组的DNA甲基化测序数据。然而，因为存在多种测序技术以及多种DNA甲基化预处理的技术，这些高通量的数据的存储、处理和分析是目前DNA甲基化研究的一个难点和热点。目前常见的高通量DNA甲基化数据检测，处理和分析的流程如图1所示。

图1 高通量DNA甲基化测序数据的检测，处理和分析的方法及软件Fig.1 Methods of detection and software packages of analysis for high-throughput sequencing of DNA methylation

3.1 DNA甲基化序列数据处理的挑战

3.1.1 甲基化预处理方法的差异和测序技术的差异

MeDIP-seq和 MBD-seq只能检测某个区域的甲基化状态，而BS-Seq、RRBS方法能够测得单碱基水平的甲基化状态。不同的DNA甲基化检测方法测得的数据也存在差异，需要不同的处理和分析方法。

3.1.2 MBD-Seq、MeDIP-Seq 数据处理的挑战

MBD-Seq和MeDIP-Seq测得的序列数据可以使用Bowtie、SOAP等短序列比对软件直接比对到参考基因组上，用映射到某个区域的reads数目来反应这个区域的甲基化程度［10-11］。然而，这两种测序方法检测的区域偏向CpG密集的甲基化区域。当某个甲基化区域的CpG分散时，很有可能被视为非甲基化区域。基因组的不同区域上CpG密度分布是不均匀的，因而需要开发新的生物信息学方法来校正，以获取基因组范围内准确的甲基化水平。

3.1.3 BS-Seq、RRBS 数据处理的挑战

BS-Seq和RRBS可以直接测得单个胞嘧啶的甲基化状态，准确性很高。然而，因为经过重亚硫酸盐转换之后，DNA的序列发生了改变(C变成了T，mC和其他碱基保持不变)，不能够直接比对到参考基因组上。另外，与Illumina直观的碱基序列不同，SOLiD测序将reads利用颜色空间进行编码，将每一个碱基与它邻近的碱基用一种颜色表示。碱基序列比对的工具不适用于SOLiD测序产生的序列。

3.2 DNA甲基化序列数据处理分析的研究现状

研究者已经开发的峰度探测软件包括MACS，USeq，PeakSeq，FindPeaks，BayesPeak 等，其中 MACS是目前最常用的峰值探测工具。然而，目前仍没有专门处理MBD-seq数据的工具或软件来降低或去除CpG密度对MBD-seq产生数据的影响。

研究者基于短序列匹配算法(Bowtie，SOAP等)开发了10多种专门处理重亚硫酸盐转换后的reads的比对工具和算法，比如 Bismark，MethylCoder，BRAT，BSMAP，BS Seeker，B-SOLADA，SOCS-B，BatMeth，RMAP-BS，FadE 等［12-14］。其中，Bismark是最常用的碱基序列比对工具，FadE，BSOLADA，SOCS-B，BatMeth是可以处理颜色空间编码的reads。如表1所示。

表1 2011～2012年BS－Seq分析软件包比较Table1 Comparison of software packages for BS－Seq analysis from 2011 to 2012

4 BS-Seq的数据处理及分析

4.1 BS-Seq 的原理

BS-Seq先利用重亚硫酸盐转换将普通的胞嘧啶变为U，而甲基化的胞嘧啶保持不变，然后使用PCR扩增使得U变成T。对转换和扩增后的DNA序列进行测序，将得到的DNA序列与参考基因组进行比较。认为C-C配对(参考基因组上在某个位置上是C，测得的reads在该位置上也是C)的就是甲基化的胞嘧啶，C-T配对的是非甲基化的胞嘧啶。如图2所示。

图2 BS－Seq原理Fig.2 BS －Seq protocol

4.2 BS-Seq数据处理流程

使用BS-Seq测得的序列数据通常为fastq或fasta格式。从序列数据中获得单个胞嘧啶的甲基化水平一般包括以下几个步骤，如图3所示:

图3 BS－Seq数据处理流程Fig.3 Recommended workflow for the analysis of BS－Seq data

(1)序列的质量控制。对于真实的数据，当reads的长度增加时，测序的错误率倾向于升高。另外，reads上包含的引物会降低匹配到基因组上的准确率。因此，有时候会对序列数据进行碱基质量分数控制、修剪引物等处理。

(2)序列比对。BS-Seq产生的序列与基因组上的原始序列存在差异(普通C变为T，互补链上的G变成了A)，需要使用BS-Seq特有的序列比对软件(Bismark等)，将BS-Seq产生的序列数据比对到参考基因组上。

(3)产生甲基化水平。从reads的基因组位置中获得每个胞嘧啶的甲基化reads数和非甲基化reads数。然后使用公式M/(U+M)计算某个胞嘧啶的甲基化水平，U和M分别是在这个胞嘧啶上的非甲基化reads数和甲基化reads数。

4.3 甲基化水平的后续分析

将单个胞嘧啶上的测序信息转换成了［0，1］的DNA甲基化水平后，研究者开发了一系列的DNA甲基化数据分析工具，实现从DNA甲基化水平中寻找甲基化模式和统计学分析等功能，以方便实验生物学家进行进一步的DNA甲基化调控机制的研究。

张岩教授课题组于2012年开发了一个可视化工具CpG_MPs，可以从标准化后的DNA甲基化水平中筛选甲基化区域和非甲基化区域［15］。Altuna等人也于2012年开发了一个R包，实现了对DNA甲基化水平的样本质量可视化、差异甲基化分析、功能注释等功能［16］。

5 总结

基于二代测序技术的DNA甲基化检测方法极大地推动了DNA甲基化的研究。研究者基于这些技术产生的高通量数据开发了一系列的生物信息学工具，然而，仍然有许多问题需要解决。目前已经开发了许多种工具可以处理和分析BS-Seq数据，然而对于MBD-Seq和MeDIP-Seq，虽然也有一些工具，但却还无法解决CpG密度偏性的问题。对于BS-Seq的数据，颜色空间编码的碱基序列比对的精度和效率依然是一项挑战。

References)

［1］ LAIRD P W.Principles and challenges of genomewide DNA methylation analysis［J］.Nature reviews.Genetics，2010 ，11，191-203.

［2］ BIRD A.DNA methylation patterns and epigenetic memory［J］.Genes＆ development，2002 ，16，6-21.

［3］ GORE A，LI Z，FUNG H L，et al.Somatic coding mutations in human induced pluripotent stem cells［J］.Nature，2011，471，63-67.

［4］ SU Jianzhong，ZHANG Yan，LÜ Jie，et al.CpG_MI:a novel approach for identifying functional CpG islands in mammalian genomes［J］.Nucleic Acids Res，2009，38，e6.

［5］ ZILBERMAN D，HENIKOFF S.Genome-wide analysis of DNA methylation patterns［J］.Development，2007，134，3959-3965.

［6］ LISTER R，PELIZZOLA M，DOWEN R H，et al.Human DNA methylomes at base resolution show widespread epigenomic differences［J］.Nature，2009，462，315-22.

［7］ LISTER R，O＇MALLEY R C，TONTI-FILIPPINI J，et al.Highly integrated single-base resolution maps of the epigenome in Arabidopsis［J］.Cell，2008，133，523-36.

［8］ STADLER M B，MURR R，BURGER L，et al.DNA-binding factors shape the mouse methylome at distal regulatory regions［J］.Nature，2011，484，550.

［9］ LÜ Jie，LIU Hongbo，SU Jianzhong，et al.DiseaseMeth:a human disease methylation database［J］.Nucleic Acids Res，2012，40，D1030-1035.

［10］ LI Ruiqiang，YU Chang，LI Yingrui，et al.SOAP2:an improved ultrafast tool for short read alignment［J］.Bioinformatics，2009，25，1966-1967.

［11］ LANGMEAD B，TRAPNELL C，POP M，et al.Ultrafast and memory-efficient alignment of short DNA sequences to the human genome［J］.Genome biology，2009，10，R25.

［12］ KRUEGER F，KRECK B，FRANKE A，et al.DNA methylome analysis using short bisulfite sequencing data［J］，Nat Methods ，2012，9，145-151.

［13］ LIM J Q，TENNAKOON C，LI G，et al.BatMeth:improved mapper for bisulfate sequencing reads on DNA methylation［J］，Genome Biology，2012，13:R82.

［14］ SOUAIAIA T，ZHANG Z， CHEN T.FadE:whole genome methylation analysisformultiplesequencing platforms［J］.Nucleic Acids Res，2012，41，e14.

［15］ SU Jianzhong，YAN Haidan，WEI Yanjun，et al.CpG_MPs: identification ofCpG methylation patternsof genomic regions from high- throughput bisulfite sequencing data［J］.Nucleic Acids Res，2012，41，e4.

［16］ AKALIN A，KORMAKSSON M，LI S，et al.methylKit:a comprehensive R package for the analysis of genomewide DNA methylation profiles［J］，Genome Biology，2012，13:R87.