DNA 甲基化数据分析流程（bismark + methylKit）

DNA甲基化在维持细胞正常功能、传递基因组印记，胚胎发育、肿瘤发生等方面发挥重要作用，目前已经成为表观遗传学和表观基因组学的研究热点。针对不同的研究目的，基于高通量测序的甲基化研究方法也是日新月异，色彩缤纷，本文简单总结一下比较常规的全基因组DNA甲基化测序数据的分析流程。基本流程分为两个步骤，mapping（bismark软件）+ 全基因组水平差异甲基化分析（methylKit软件）。

 Mapping

bismark是目前WGBS或RRBS产生的甲基化数据最受欢迎的软件之一，操作简单，功能强大，bismark主要由3个步骤组成。

    a)bismark_genome_preparation. 使用Bisulfite处理样品DNA是进行DNA甲基化分析的一种常用方法，处理后会使未发生甲基化的胞嘧啶(Cytosine)转换为胸腺嘧啶(Thymine)，导致按照常规的比对方法无法将测序得到的reads比对到参考基因组上。软件Bismark首先会对参考基因组进行转换，使甲基化数据与参考基因组具有可比性。

Usage，bismark_genome_preparation [options] <path_to_genome_folder>

    b). bismark（alignment 步骤），这一步会对测序数据（每一条reads），做与参考基因组同样的处理，然后把reads向基因组上mapping，这一步会得到比对结果bam/sam文件

  bismark [options] <genome_folder> {-1 <mates1> -2 <mates2> |<singles>}

    c). Bismark bismark_methylation_extractor,，这是bismark的最后一步，根据bam/sam文件提取每条reads中每个C位点的甲基化状态。结果会根据C位点所在的序列context分别存在CpG_context_bismark.txt，CHG_context_bismark.txt，CHH_context_bismark.txt这三个文件中。同时根据CpG_context_bismark.txt报告的CpG甲基化信息，bismark会计算出全基因组内每个CpG位点的甲基化水平，存放在bismark.CpG_report.txt文件中。

USAGE: bismark_methylation_extractor [options] <filenames>

运行完上述3个步骤，我们就获得了最为基本的DNA甲基化水平的信息，可以以此为基础，进行下游的分析，例如，对比不同样本之间的差异甲基化，等等。

2.差异甲基化分析（methylKit 软件）

methylKit是基于R语言开发的程序包，对于熟悉R语言的小伙伴来说，同样是简单粗暴，功能强大。（此处略过其安装过程，可参考https://github.com/al2na/methylKit）

    a). 数据读入，methylKit可以利用其中的read.bismark函数直接读取bismark产生的测序数据（用法如图所示），不过并不建议这么做，因为R语言运行速度较慢，读取sam文件会非常耗时。我们可以直接根据bismark产生的bismark.CpG_report.txt文件，手动写一个脚本，把数据处理成methylKit需要的格式（如图）

    然后利用read函数读取文件内容（如图）

b)绘制数据的基本统计图

c)所有样本合并成同一个数据集

d)统计样本之间的相关性

e)样本聚类分析

f)PCA分析

g) 差异甲基化分析（可以以位点为单位或者以一定长度的窗口为单位）

以位点为单位:

以窗口为单位，首先要整理出每个样本对应同一个窗口的DNA甲基化水平，可以自己定义步长和窗口大小。

然后一样的过程计算差异甲基化的窗口

此时要注意把参数meth 换为tiles。

h)对差异甲基化区域进行注释

methylKit包可以帮助我们完成所有的关于DNA甲基化数据的基本分析，是不是很强大。而下游更具体的分析，因为每个人的关注点不同，因此没有一个统一的流程，需要根据自己的需求自行设计方案。

欢迎关注