高通量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.