Whole Genome DNA Methylation Ana

一、介绍

首先声明，本文注重流程的梳理和生物学意义的解答，具体软件的原理和参数的运用请参照MethPipe软件手册（和我未来的翻译）

本文使用Brachypodium distachyon作为模式植物

二、本文所需的软硬件

该部分后续于代码中将直接安装，无需于此步骤手动处理

三、正式开始甲基化测序数据分析

3.1 生成基因组索引文件（3.1和3.2步骤使用的软件名为WALT）

首先是把B. distachyon Bd21 基因组和chloroplast（叶绿体）测序文件融合成一个文件

$zcatBdistachyon_314_v3.0.fa.gz|cat-chloroplast.fa> Bdistachyon_314_v3.0_ch.fa

然后用makedb函数生成上述文件的索引

$ makedb -c Bdistachyon_314_v3.0_ch.fa -o Bdistachyon_314_v3.0_ch.dbindex

3.2 将BS-seq的测序片段比对至已建立好索引的参考基因组上

Mapping

注意：

黄色高亮：在双端测序中，Read_1文件（即5'端到3'端）碱基T富集，而Read_2文件（即3'端到5'端）碱基A富集。

红色下划线：walt函数的参数-c用于修建Illumina标准接头（测序时人为添加的）

最后一段to-mr函数：walt函数可输出.mr文件或.sam文件，此函数用于转换

3.3 去除重复序列（3.3、3.4、3.5步骤使用的软件名为MethPipe）

生物学原理：如果测序片段之间具有相同的序列并比对至基因组相同的位置，那么这个现象很有可能是PCR扩增导致的，因此需要在后续差异分析之前去除这些重复。

在去除重复之前，需对上步得到的MR文件进行排序（染色体、起点、止点、链），这一步很重要，很多其他的文件类型如.bed等也对文件内的数据格式有要求。

Sort和Duplicate-remover

然后对该文件去除重复，其后的参数请详见MethPipe软件（该软件的手册已附于文末）

3.4 估计重亚硫酸盐转换的比例（即该实验处理的完全性）

生物学原理：叶绿体基因组被用于作为对照组（control），因为我们共同认为叶绿体基因组中的胞嘧啶（C）均未被甲基化。

首先用grep函数，从上步中已比对好的测序数据中，将叶绿体基因组的测序数据单独拎出。

Grep

然后用bsrate函数，估计重亚硫酸盐反应的转换率，越靠近1表示转换的越完全。

Bsrate

3.5 计算甲基化水平和其他相关的统计数据——统计单个碱基的甲基化水平

函数methcounts：对样本中所有的胞嘧啶（C）统计其甲基化水平。

Methcounts

函数levles：将上述统计结果进行统计学分析，该函数主要计算如下图所示。

Levels

到此为止，甲基化测序数据的上游处理已经结束，该文章同时给出了三种R包进行下游处理，分别为methylKit、EnrichedHeatmap和methylPipe，这三个包加上上游处理需要用到的MethPipe软件和WALT软件应该都可以从Github上下载。

最后提一句的是，用macOS的RStudio中从Github上下载软件需要安装Command Line插件（该插件应该已经镶嵌于XCODE中）但如果没有，请如下图操作。

Command Line