DNA甲基化
2019-07-23 本文已影响35人
六六_ryx
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)
甲基化芯片的分析可以参考：https://www.jianshu.com/p/6411e8acfab3
