一、介绍
首先声明，本文注重流程的梳理和生物学意义的解答，具体软件的原理和参数的运用请参照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
网友评论