DNA甲基化

DNA甲基化是表观遗传学研究中的一个重要一个分支，符合表观遗传学DNA的序列信息没有改变的特地点，而是DNA的碱基发生了化学修饰。最常见的一种修饰如5-甲基胞嘧啶（5-methylcytosine, 5mC）。相关研究有关注于甲基化与基因表达调控、基因印记、X染色体失活、胚胎发育、肿瘤发生等。

甲基化的研究技术

总体上来说可以分为两类：一种是测序，一种是芯片。

测序
如全基因组DNA甲基化测序（Whole Genome Bisulfite Sequencing，WGBS）和简化甲基化测序 (Reduced representation bisulfite sequencing, RRBS)，前者是利用 Bisulfite（重亚硫酸盐）将未甲基化的C碱基转换为U，而甲基化的C碱基不变的原理，从而对甲基化信号检测。后者是根据酶切 (Msp I) 富集启动子及CpG岛区域，并进行Bisulfite测序，同时实现DNA甲基化状态检测的高分辨率和测序数据的高利用率。

BS-seq原理
甲基化芯片
Illumina的Infinium BeadChip芯片，包括HumanMethyation450（450K）和MethylationEPIC（850K）

甲基化分析流程

上游分析的主要目的提取甲基化位点，得到差异甲基化位点和甲基化区域（DML/DMR, DNA methylation loci/region）,下游分析比较灵活，甲基化分析通常与其他组学联合分析，如：

联合转录组，对DML/DMR区域内的基因进行功能富集分析（可用的软件有TopGO ，GSEABase ， Enrichr），以及与非编码RNA的整合分析等
联合组蛋白修饰，一是借助组蛋白修饰可以注释DML/DMR基因组分布，是否在启动子或增强子等区域
联合基因组，分析与甲基化与突变、变异的关系

STEP1:比对，提取甲基化位点

对于BS-seq甲基化位点的提取，常用的软件有Bismark，bismark可以将bisulfite处理的reads比对到参考基因组上，同时提取甲基化位点。输出结果可以直接导入到基因组浏览器（如IGV, SeqMonk）中查看甲基化水平。主要特点可以概括为以下几个方面：

一步实现Bisulfite 比对和提取甲基化位点和区域
支持单端和双端测序
支持ungapped, gapped or spliced alignments
输出区分胞嘧啶甲基化在CpG, CHG和CHH
流程代码参考：https://github.com/FelixKrueger/Bismark/tree/master/Docs

使用方法
主要分为3步：

1 . 准备参考基因组，需要将参考基因组进行bisulfite 转换和建索引, 建索引是基于bowtie2和hisat2 bowtie2-build or hisat2-build
USAGE: bismark_genome_preparation [options] <path_to_genome_folder>

bismark_genome_preparation --path_to_aligner /usr/bin/bowtie2/ --verbose /data/genomes/homo_sapiens/GRCh38/

2 . 比对
比对前需要先进行质控，默认比对模式是bowtie2，输出文件格式默认是bam文件，标准的比对方式设置multi-seed 的长度为20bp，即-L 20，0 错配：-N 0。
USAGE:bismark [options] --genome <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}
单端测序：

bismark --bowtie2 -N 0 -L 20  -o mapping-out ref data.fastq
# ref:是参考基因组的文件夹
# -o: mapping-out输出文件夹
# --bowtie2：默认模式，可选参数

3 .提取甲基化位点
甲基化输出文件中的符号的含义：点.和小写的z，x,h,u都表示无甲基化，代表的类型分别是CpG，CHG，CHH，CN or CHN，大写的z，x,h,u都表示有甲基化。

在提取甲基化位点前需要先去重, 不过对RRBS, amplicon以及靶向测序不建议去重。
USAGE: ./deduplicate_bismark [options] filename(s)

deduplicate_bismark -s mapping-out/test_data_bismark_bt2.bam

call DML/DMR

mkdir CpGout
bismark_methylation_extractor -s --gzip --bedGraph --buffer_size 1G \
--cytosine_report --genome_folder ref/ test_data_bismark_bt2.deduplicated.bam

输出结果：

第一列测序信息
第二列甲基化状态，+代表甲基化，- 代表未甲基化
第三列代表chromosome
第四列代表location
第五列代表methylation call

STEP2：差异甲基化区域分析

差异甲基化位点和区域（DML/DMR）的区别

差异甲基化区域的分析工具很多，这里介绍两个R包，一个是edmr, 一个是DSS。

edmr

edmr是基于双峰正态分布模型和区域甲基化分析的成本函数优化DMR分析https://github.com/ShengLi/edmr](https://github.com/ShengLi/edmr。

# Step 1. Load add-on packages and example data

library(edmr)
library(GenomicRanges)
library(IRanges)
library(mixtools)
library(data.table)
data(edmr)
# Step 2. myDiff evalution and plotting

# fitting the bimodal normal distribution to CpGs distribution
myMixmdl=myDiff.to.mixmdl(myDiff, plot=T, main="example")

# plot cost function and the determined distance cutoff
plotCost(myMixmdl, main="cost function")

# Step 3. Calculate DMRs

# calculate all DMRs candidate
mydmr=edmr(myDiff, mode=1, ACF=TRUE)

# further filtering the DMRs
mysigdmr=filter.dmr(mydmr)

## annotation
# get genebody annotation GRangesList object
#genebody=genebody.anno(file="http://edmr.googlecode.com/files/hg19_refseq_all_types.bed")
genebody.file=system.file("extdata", "chr22.hg19_refseq_all_types.bed.gz", package = "edmr")
genebody=genebody.anno(file=genebody.file)

# plot the eDMR genebody annotation
plotdmrdistr(mysigdmr, genebody)

# get CpG islands and shores annotation
#cpgi=cpgi.anno(file="http://edmr.googlecode.com/files/hg19_cpgisland_all.bed")
cpgi.file=system.file("extdata", "chr22.hg19_cpgisland_all.bed.gz", package = "edmr")
cpgi=cpgi.anno(file=cpgi.file)

# plot the eDMR CpG islands and shores annotation
plotdmrdistr(mysigdmr, cpgi)

# prepare genes for pathway analysis with significant DMRs at its promoter regions 
dmr.genes=get.dmr.genes(myDMR=mysigdmr, subject=genebody$promoter, id.type="gene.symbol")
dmr.genes

DSS

DSS是基于贝塔二项分布，可以用于分析DML/DMRs)，包括三个模块： -两组比较（有重复），两组比较（无重复），多组比较，参考https://www.bioconductor.org/packages/release/bioc/html/DSS.html。
DSS的输入格式如下：

第一列是染色体号
第二列是位置
第三列总的read counts
第四列是甲基化的read counts
第三列的数字要大于或等于第四列，要保证第三列和第四列数据类型是数字格式

library(DSS)
# 创建BSseq对象
if(T){
 sn[c(1,11,2,12)]
 BSobj <- makeBSseqData(allDat[c(1,11,2,12)],sn[c(1,11,2,12)])[1:5000,]
 BSobj 
 save(BSobj,file = 'group-BSobj.Rdata')
 dmlTest <- DMLtest(BSobj,   group1=c("A0R_d0_rep1","A0R_d0_rep2"),
 group2=c("A3R_d0_rep1","A3R_d0_rep2"),smoothing=F)
 head(dmlTest) 
}
##  call DML/DMR
# using callDML function to call DML
# call DML
dmls = callDML(dmlTest, p.threshold=0.001)
head(dmls)
# call DMR
dmrs = callDMR(dmlTest, p.threshold=0.01)
head(dmrs)
# visulization
showOneDMR(dmrs[1,], BSobj)

甲基化芯片的分析可以参考：//www.greatytc.com/p/6411e8acfab3