高通量DNA甲基化數(shù)據(jù)的處理和分析方法

王心宇，許穎出，劉洪波，王 芳，張 巖，蘇建忠

(哈爾濱醫(yī)科大學(xué)生物信息科學(xué)與技術(shù)學(xué)院，黑龍江哈爾濱150080)

DNA甲基化是重要的表觀遺傳學(xué)修飾之一，以往的研究表明，DNA甲基化在細胞發(fā)育和分化、調(diào)控基因表達、X染色體失活、基因沉默、疾病的發(fā)生等方面扮演著重要的角色［1-3］。在真核生物中，通常是CpG二核苷酸中胞嘧啶的第五個碳原子上發(fā)生了甲基化(5mC)，胞嘧啶甲基化也可能會發(fā)生在CHG和CHH(H是除G外的任意一種核苷酸)上。全基因組的甲基化水平呈現(xiàn)雙峰分布，而且低甲基化的區(qū)域多數(shù)是在CpG二核苷酸聚集區(qū)域(CpG島)［4］。以往的研究發(fā)現(xiàn)，位于啟動子區(qū)域的高甲基化的CpG島與基因的沉默有關(guān)，可能是因為DNA甲基化阻礙轉(zhuǎn)錄因子結(jié)合而直接抑制了基因轉(zhuǎn)錄。在過去的幾十年里，由于實驗技術(shù)和費用的限制，DNA甲基化的數(shù)據(jù)往往只檢測了基因組的局部區(qū)域，而且是低通量的數(shù)據(jù)。

二代測序技術(shù)的發(fā)展極大地推動了表觀遺傳調(diào)控機制的研究，基于二代測序技術(shù)發(fā)展起來的DNA甲基化的檢測技術(shù)為DNA甲基化的研究提供了大量的高通量、全基因組的DNA甲基化數(shù)據(jù)。這些高通量數(shù)據(jù)的產(chǎn)生使得DNA甲基化研究的重點由目標基因DNA甲基化的檢測轉(zhuǎn)移到了全基因組DNA甲基化高通量數(shù)據(jù)的檢測、存儲、處理和分析上。近幾年，研究者構(gòu)建了多個DNA甲基化數(shù)據(jù)庫，開發(fā)了大量的DNA甲基化高通量數(shù)據(jù)的處理和分析工具，使得深入的表觀遺傳調(diào)控機制的研究成為可能。

1 基于二代測序技術(shù)的DNA甲基化檢測技術(shù)

1.1 DNA甲基化預(yù)處理方法

甲基化后的胞嘧啶(5 mC)與普通的胞嘧啶(C)在DNA序列上并無差異，如果直接使用DNA測序，將無法區(qū)分測得的胞嘧啶C是C還是5 mC。所以檢測DNA甲基化需要首先對待檢測的DNA序列中胞嘧啶進行預(yù)處理，將非甲基化的胞嘧啶C與甲基化的胞嘧啶5 mC區(qū)分開來，目前的DNA甲基化預(yù)處理方式主要分為三種:

(1)限制性內(nèi)切酶法(Endonuclease digestion)

限制性內(nèi)切酶法是指利用甲基化限制性內(nèi)切酶(HpaII，MspI和HhaI等)在各自的識別位點對甲基化的胞嘧啶有不同的敏感性來檢測CpG的甲基化［5］。限制性內(nèi)切酶法結(jié)合二代測序的技術(shù)有MRE-seq，MCA-seq，MSCC 和 HELP-seq。盡管限制性內(nèi)切酶測序法成本低、高效，然而由于檢測的CpG位點局限于酶切位點附近，基因組覆蓋率低，另外還存在CpG偏好性、酶切不完全導(dǎo)致的假陽性等問題，使用這種方法檢測DNA甲基化的研究越來越少。

(2)親和純化法(Affinity enrichment)

親和純化是利用甲基化CpG結(jié)合蛋白(MBD)或者對5mC特異的抗體來親和提純甲基化區(qū)域。MeDIP-seq和 MBD-seq是最常用的兩種結(jié)合親和純化和二代測序技術(shù)的DNA甲基化檢測方法?；跍y序的親和純化法能夠快速、低成本地檢測全基因組范圍內(nèi)的甲基化水平，然而它只能獲得區(qū)域的甲基化水平，特別是MeDIP-seq偏向于CpG富集的區(qū)域，分散的低密度的甲基化位點可能被識別成非甲基化區(qū)域，目前還沒有能夠去除掉這種偏性的生物信息學(xué)方法。

(3)重亞硫酸鹽轉(zhuǎn)換法(Bisulphite conversion)

重亞硫酸鹽轉(zhuǎn)換結(jié)合二代測序技術(shù)是目前最精準的DNA甲基化檢測方法，能夠檢測單堿基水平的甲基化狀態(tài)，被稱為DNA甲基化檢測的“金標準”。對基因組中未發(fā)生甲基化的胞嘧啶進行重亞硫酸鹽處理，將其轉(zhuǎn)換成U，經(jīng)PCR擴增后變成T，重亞硫酸鹽轉(zhuǎn)換對甲基化的胞嘧啶不起作用。通過結(jié)合二代測序，即可繪制出單堿基分辨率的全基因組DNA甲基化圖譜。目前常用的重亞硫酸鹽轉(zhuǎn)換結(jié)合二代測序技術(shù)的DNA甲基化檢測技術(shù)有BS-seq和RRBS等。

1.2 二代測序技術(shù)

在使用DNA甲基化預(yù)處理區(qū)分出未甲基化的胞嘧啶和甲基化的胞嘧啶后，再使用二代測序技術(shù)檢測DNA序列，來獲取胞嘧啶上的甲基化狀態(tài)。

目前二代測序技術(shù)主要分為三個平臺:Roche、Illumina、SOLiD。其中每種測序平臺又擁有多種系統(tǒng)，比如Illumina就有HiSeq、GAIIx等系統(tǒng)。不同的測序技術(shù)在測得的read長度、精確性、通量都有差異，適用于不同的研究目的需要。

1.3 高通量DNA甲基化檢測技術(shù)進展

結(jié)合二代測序技術(shù)和DNA甲基化預(yù)處理的DNA甲基化檢測方法，在近幾年獲得了大量的全基因組的DNA甲基化測序數(shù)據(jù)。

國外很多實驗室產(chǎn)生了大量、精準的高通量DNA甲基化數(shù)據(jù)，例如，Lister等人于2008年檢測的擬南芥全基因組甲基化譜和2009年測得的人類全基因組甲基化譜［6-7］，Stadler等人于2011年測定了小鼠胚胎干細胞和神經(jīng)前體細胞的全基因組甲基化譜等［8］。國內(nèi)近年來也產(chǎn)生了大量的高通量DNA甲基化數(shù)據(jù)，例如，2010年，中科院昆明研究所，華大基因和上海交通大學(xué)癌癥表觀遺傳中心等九家科研機構(gòu)聯(lián)合測定了桑蠶的單堿基水平的DNA甲基化譜，王俊教授課題組測定的人類完全分化的血細胞的全基因組DNA甲基化譜等。這些全基因組水平的DNA甲基化數(shù)據(jù)為表觀遺傳調(diào)控機制的研究提供了數(shù)據(jù)資源。

2 DNA甲基化數(shù)據(jù)儲存和可視化

目前研究者構(gòu)建了各種各樣的數(shù)據(jù)庫來存儲世界范圍的各大實驗室和科研機構(gòu)產(chǎn)生的高通量DNA甲基化數(shù)據(jù)，便于數(shù)據(jù)的查詢、下載、可視化分析及全球化的資源共享。從第一個DNA甲基化的公共數(shù)據(jù)庫MethDB由Grunau等人于2001年構(gòu)建以來，已有多個和DNA甲基化相關(guān)的數(shù)據(jù)庫被開發(fā)，例如，NCBI的存儲表觀遺傳修飾數(shù)據(jù)的Epigenomics，主要包括DNA甲基化、組蛋白修飾和非編碼RNA等數(shù)據(jù)。PubMeth是結(jié)合文本的基因注釋信息的DNA甲基化數(shù)據(jù)庫。DiseaseMeth儲存72種人類疾病相關(guān)的DNA甲基化的數(shù)據(jù)庫，并實現(xiàn)了統(tǒng)計學(xué)分析及可視化［9］。

3 高通量DNA甲基化數(shù)據(jù)的處理和分析

結(jié)合二代測序技術(shù)和DNA甲基化預(yù)處理的方法，在近幾年產(chǎn)生了大量的全基因組的DNA甲基化測序數(shù)據(jù)。然而，因為存在多種測序技術(shù)以及多種DNA甲基化預(yù)處理的技術(shù)，這些高通量的數(shù)據(jù)的存儲、處理和分析是目前DNA甲基化研究的一個難點和熱點。目前常見的高通量DNA甲基化數(shù)據(jù)檢測，處理和分析的流程如圖1所示。

圖1 高通量DNA甲基化測序數(shù)據(jù)的檢測，處理和分析的方法及軟件Fig.1 Methods of detection and software packages of analysis for high-throughput sequencing of DNA methylation

3.1 DNA甲基化序列數(shù)據(jù)處理的挑戰(zhàn)

3.1.1 甲基化預(yù)處理方法的差異和測序技術(shù)的差異

MeDIP-seq和 MBD-seq只能檢測某個區(qū)域的甲基化狀態(tài)，而BS-Seq、RRBS方法能夠測得單堿基水平的甲基化狀態(tài)。不同的DNA甲基化檢測方法測得的數(shù)據(jù)也存在差異，需要不同的處理和分析方法。

3.1.2 MBD-Seq、MeDIP-Seq 數(shù)據(jù)處理的挑戰(zhàn)

MBD-Seq和MeDIP-Seq測得的序列數(shù)據(jù)可以使用Bowtie、SOAP等短序列比對軟件直接比對到參考基因組上，用映射到某個區(qū)域的reads數(shù)目來反應(yīng)這個區(qū)域的甲基化程度［10-11］。然而，這兩種測序方法檢測的區(qū)域偏向CpG密集的甲基化區(qū)域。當某個甲基化區(qū)域的CpG分散時，很有可能被視為非甲基化區(qū)域?；蚪M的不同區(qū)域上CpG密度分布是不均勻的，因而需要開發(fā)新的生物信息學(xué)方法來校正，以獲取基因組范圍內(nèi)準確的甲基化水平。

3.1.3 BS-Seq、RRBS 數(shù)據(jù)處理的挑戰(zhàn)

BS-Seq和RRBS可以直接測得單個胞嘧啶的甲基化狀態(tài)，準確性很高。然而，因為經(jīng)過重亞硫酸鹽轉(zhuǎn)換之后，DNA的序列發(fā)生了改變(C變成了T，mC和其他堿基保持不變)，不能夠直接比對到參考基因組上。另外，與Illumina直觀的堿基序列不同，SOLiD測序?qū)eads利用顏色空間進行編碼，將每一個堿基與它鄰近的堿基用一種顏色表示。堿基序列比對的工具不適用于SOLiD測序產(chǎn)生的序列。

3.2 DNA甲基化序列數(shù)據(jù)處理分析的研究現(xiàn)狀

研究者已經(jīng)開發(fā)的峰度探測軟件包括MACS，USeq，PeakSeq，F(xiàn)indPeaks，BayesPeak 等，其中 MACS是目前最常用的峰值探測工具。然而，目前仍沒有專門處理MBD-seq數(shù)據(jù)的工具或軟件來降低或去除CpG密度對MBD-seq產(chǎn)生數(shù)據(jù)的影響。

研究者基于短序列匹配算法(Bowtie，SOAP等)開發(fā)了10多種專門處理重亞硫酸鹽轉(zhuǎn)換后的reads的比對工具和算法，比如 Bismark，MethylCoder，BRAT，BSMAP，BS Seeker，B-SOLADA，SOCS-B，BatMeth，RMAP-BS，F(xiàn)adE 等［12-14］。其中，Bismark是最常用的堿基序列比對工具，F(xiàn)adE，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的數(shù)據(jù)處理及分析

4.1 BS-Seq 的原理

BS-Seq先利用重亞硫酸鹽轉(zhuǎn)換將普通的胞嘧啶變?yōu)閁，而甲基化的胞嘧啶保持不變，然后使用PCR擴增使得U變成T。對轉(zhuǎn)換和擴增后的DNA序列進行測序，將得到的DNA序列與參考基因組進行比較。認為C-C配對(參考基因組上在某個位置上是C，測得的reads在該位置上也是C)的就是甲基化的胞嘧啶，C-T配對的是非甲基化的胞嘧啶。如圖2所示。

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

4.2 BS-Seq數(shù)據(jù)處理流程

使用BS-Seq測得的序列數(shù)據(jù)通常為fastq或fasta格式。從序列數(shù)據(jù)中獲得單個胞嘧啶的甲基化水平一般包括以下幾個步驟，如圖3所示:

圖3 BS－Seq數(shù)據(jù)處理流程Fig.3 Recommended workflow for the analysis of BS－Seq data

(1)序列的質(zhì)量控制。對于真實的數(shù)據(jù)，當reads的長度增加時，測序的錯誤率傾向于升高。另外，reads上包含的引物會降低匹配到基因組上的準確率。因此，有時候會對序列數(shù)據(jù)進行堿基質(zhì)量分數(shù)控制、修剪引物等處理。

(2)序列比對。BS-Seq產(chǎn)生的序列與基因組上的原始序列存在差異(普通C變?yōu)門，互補鏈上的G變成了A)，需要使用BS-Seq特有的序列比對軟件(Bismark等)，將BS-Seq產(chǎn)生的序列數(shù)據(jù)比對到參考基因組上。

(3)產(chǎn)生甲基化水平。從reads的基因組位置中獲得每個胞嘧啶的甲基化reads數(shù)和非甲基化reads數(shù)。然后使用公式M/(U+M)計算某個胞嘧啶的甲基化水平，U和M分別是在這個胞嘧啶上的非甲基化reads數(shù)和甲基化reads數(shù)。

4.3 甲基化水平的后續(xù)分析

將單個胞嘧啶上的測序信息轉(zhuǎn)換成了［0，1］的DNA甲基化水平后，研究者開發(fā)了一系列的DNA甲基化數(shù)據(jù)分析工具，實現(xiàn)從DNA甲基化水平中尋找甲基化模式和統(tǒng)計學(xué)分析等功能，以方便實驗生物學(xué)家進行進一步的DNA甲基化調(diào)控機制的研究。

張巖教授課題組于2012年開發(fā)了一個可視化工具CpG_MPs，可以從標準化后的DNA甲基化水平中篩選甲基化區(qū)域和非甲基化區(qū)域［15］。Altuna等人也于2012年開發(fā)了一個R包，實現(xiàn)了對DNA甲基化水平的樣本質(zhì)量可視化、差異甲基化分析、功能注釋等功能［16］。

5 總結(jié)

基于二代測序技術(shù)的DNA甲基化檢測方法極大地推動了DNA甲基化的研究。研究者基于這些技術(shù)產(chǎn)生的高通量數(shù)據(jù)開發(fā)了一系列的生物信息學(xué)工具，然而，仍然有許多問題需要解決。目前已經(jīng)開發(fā)了許多種工具可以處理和分析BS-Seq數(shù)據(jù)，然而對于MBD-Seq和MeDIP-Seq，雖然也有一些工具，但卻還無法解決CpG密度偏性的問題。對于BS-Seq的數(shù)據(jù)，顏色空間編碼的堿基序列比對的精度和效率依然是一項挑戰(zhàn)。

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，F(xiàn)UNG 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，F(xiàn)RANKE 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.