干货|全基因组甲基化测序分析（一）

DNA甲基化是研究最为广泛的表观遗传修饰，它被认为在许多生物学过程和疾病中有着重要的作用。随着测序技术的快速发展，二代测序与重亚硫酸盐处理基因组DNA相结合产生了可以在全基因组单碱基水平上测量DNA甲基化水平的全基因组重亚硫酸盐测序（whole-genome bisulfite sequencing，WGBS）技术。然而紧随而来的，是来自生物信息学分析上的挑战，本文主要向大家介绍分析WGBS数据的第一步——基因组拼接（mapping）。

WGBS的流程与全基因组DNA测序主要区别于DNA文库的构建。以MethlyC-seq为例，将DNA片段化后收集特定长度的片段并修复DNA末端，再将双端加上3′-dAMP然后连接上甲基化的接头，接着用重亚硫酸盐处理DNA使未甲基化的胞嘧啶（C）转化为尿嘧啶（U），最后进行低循环数的PCR扩增后完成建库。

具体流程见图1：

图1：MethylC-seq 文库构建流程。

建库完后就是上机测序，得到原始数据（fastq文件）后，完成一定的质量控制后就可以进行mapping了。由上述可知WGBS在构建测序文库时DNA受到重亚硫酸盐的处理，非甲基化的C被转化成U，并在随后的PCR扩增中转化为胸腺嘧啶（T），而甲基化的C不会被影响，通过这一原理可以确定胞嘧啶的甲基化状态。

目前主要有两种用于比对WBGS测序fastq数据的方法。第一种three-letter法是将参考基因组和测序reads里所有的C转化成T，从而只利用三种碱基进行比对。目前利用最为普遍的three-letter法比对软件有Bismark、BS-Seeker2等。例如，Bismark软件首先将参考基因组的C转化成T（对于另一条DNA链则是将G转换成A）建立索引，再将测序reads也进行同样的转换，然后利用Bowtie、Bowtie2为内核进行比对，最后根据比对的结果计算出每个胞嘧啶的甲基化水平。第二种wild-card法是将参考基因组里所有的C转化成字符Y，在比对时Y可以与C和T比对上。其中，BSMAP是最为流行的wild-card法比对软件之一，其利用SOAP软件进行比对。

可以看出three-letter法类似于比对基因组测序数据的方法，只是少了一种碱基，但这造成了序列特异性的降低，从而导致一部分reads无法比对到基因组上。而在wild-card法中，含有甲基化胞嘧啶的reads可以含有四种碱基，而含有非甲基化胞嘧啶的reads可能只含有三种碱基，从而导致序列特异性的降低，使含有非甲基化胞嘧啶的reads比对率降低，在计算甲基化水平时引入偏倚，造成甲基化水平的升高，但是其整体的比对率高于three-letter法。具体原理见图2。

图2：WGBS数据比对的两种方法比较。

下面我们和大家分享Bismark快速上手教程：

首先需要将感兴趣的基因组建立索引以便Bowtie2用于比对，这一步对于同一基因组只需要做一次。

使用方法：bismark_genome_preparation [options] <path_to_genome_folder>

具体例子：/data/bismark_v0.16.1/bismark_genome_preparation --path_to_bowtie /usr/local/bowtie2/ /data/genomes/hg19/

重要参数：

--bowtie1、--bowtie2（Bowtie 1 与Bowtie 2所用到的索引不兼容，在建立索引时需要特定的指定。）

path_to_genome_folder（含有参考基因组序列（如hg19.fasta）的文件夹路径。）

建立索引后Bismark将调用Bowtie进行序列比对。

使用方法：bismark [options] <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}

具体例子：/data/bismark_v0.16.1/bismark --temp_dir /data/temp_bismark/ -N 1 -L 50 /data/genomes/hg19 input.fastq

重要参数：

-I（pair-end测序数据比对时最小的有效插入片段长度。如果I被指定为70，那么两个长度为20-bp的用于比对的序列之间在基因组上最小允许相隔30-bp，若两个序列在基因组上相距小于30-bp，则这两个序列不会被比对到基因组上。）

-X（pair-end测序数据比对时最大的有效插入片段长度。如果X被指定为100，那么两个长度为20-bp的用于比对的序列之间在基因组上最大允许相隔60-bp，若两个序列在基因组上相距大于60-bp，则这两个序列不会被比对到基因组上。通常I使用默认值0即可，而X则需要根据pair-end测序建库时平均的片段插入长度决定。）

--multicore、-p（使用者可根据自己的计算机或者服务器的核数与内存并行运行bismark以提高比对效率。）

--non_directional、--pbat（这两个参数与测序建库的类型相关，使用者应清楚自己的测序建库类型。）

--temp_dir（bismark在进行比对时会产生非常大的临时文件，需要将临时文件储存在指定的空间较大的地方。）

--bowtie1（bismark默认是调用Bowtie 2进行比对，想用Bowtie 1比对需要指定。）

-N（在进行multiseed比对时一个seed被允许的最大错配数。N越大则比对的越慢（通常会慢很多），但是会使比对的灵敏度提高。如果使用Bowtie 1，类似的参数见-n。）

-L（在进行multiseed比对时一个seed的长度。L越大则比对的越慢，但是会使比对的灵敏度提高。如果使用Bowtie 1，类似的参数见-l。）

<genome_folder>（这个文件夹里包含了感兴趣的基因组和上一步建立索引时所产生的文件夹Bisulfite_Genome。）

-1、-2、<singles>（如果比对的是pair-end测序的数据，则-1与-2分别对应一对原始数据。若比对的是single-end测序数据，则输入的文件只有一个。）

比对完后对每个位点进行甲基化水平的计算。

使用方法：bismark_methylation_extractor [options] <filenames>

具体例子：/data/bismark_v0.16.1/bismark_methylation_extractor -s --no_header --multicore 10 --ignore 1 --cytosine_report --genome_folder /data/genomes/hg19/ input.bam

重要参数：

-s、-p（输入的数据是single-end测序还是pair-end测序的比对结果。）

--no_overlap（若是pair-end测序，则理论上一对序列可能存在重叠部分，重叠部分的位点会被计算两次甲基化水平，造成偏倚。启用这个参数会只计算pair 1上重叠部分的位点甲基化水平。）

--ignore、--ignore_r2、--ignore_3prime、--ignore_3prime_r2（忽略pair1（pair2）的最前（最后）几个位点的信息，可有效去除序列5’端和3’端甲基化水平的偏倚，详见说明书M-bias相关说明。）

软件许多参数使用默认值即可，使用者可根据Bismark官网（http://www.bioinformatics.babraham.ac.uk/projects/bismark/）的说明书根据自己的个性化需求修改参数。

参考文献：

Urich M A, Nery J R, Lister R, et al. MethylC-seq library preparation for base-resolution whole-genome bisulfite sequencing[J]. Nature Protocols, 2015, 10(3):475.

Bock C. Analysing and interpreting DNA methylation data[J]. Nature Reviews Genetics, 2012, 13(10):705-19.

Krueger F, Andrews S R. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications.[J]. Bioinformatics, 2011, 27(11):1571-2.