0.质控
trim_galore -q 20 --phred33 --illumina --stringency 1 -e 0.1 \
--length 36 --rrbs --cores 4 --fastqc -o ./cut_adapter \
--paired R1.fastq R2.fastq
-q 设定phred quality阈值，默认20（99%的read质量），如果测序深度较深，可以设定25
--phred33 默认记分方式，代表Q+33=ASCII码的方式来记分方式
--stringency 设定可以忍受的前后adapter重叠的碱基数，默认是1
--illumina 指定要去除的接头序列，可先前运行一下fastqc来确定接头序列类型
--length 设定长度阈值，小于此长度会被抛弃
--rrbs 指定输入文件是MspI消化的RRBS样本（识别位点：）CCGG
--fastqc 修剪完成后，以默认模式在FastQ文件上运行FastQC
-o 输出目录
-e 设定默认质量控制数，默认是0.1，即ERROR rate大于10%的read 会被舍弃，如果添加了--paired参数则会舍弃一对reads
--paired 对于双端结果，一对reads中若一个read因为质量或其他原因被抛弃，则对应的另一个read也抛弃。
<filename> #如果是采用illumina双端测序的测序文件，应该同时输入两个文件。
1.建立亚硫酸盐基因组index
bismark_genome_preparation --bowtie2 基因组所在目录（目录下应含有fa或fasta格式的文件）
2. 比对
bismark --bowtie2 -p 4 -N 0 -L 20 --quiet --un --ambiguous 亚硫酸盐基因组index所在目录 -1 fastq1 -2 fastq2 -o 输出目录
-p 线程数
-N 允许错配数目
-L 比对时建立的seed大小
--quiet 不在控制台输出比对流程信息
--un 输出不匹配序列信息到_unmapped_reads.fq.gz中
--ambiguous 输出多处匹配序列信息到_ambiguous_reads.fq中，不指定则输出到_unmapped_reads.fq.gz中
-o 输出目录
--sam 指定输出的格式为sam，无此参数则输出bam格式文件
输出结果有两个重要文件，一个bismark_bt2_PE_report.txt，一个bismark_bt2_pe.bam
3-1.去重（可选，建议对WGBS使用此步骤，但不应将其用于代表性较低的文库，例如RRBS，扩增子或靶标富集文库）
deduplicate_bismark -p --bam 输入的bam文件 -o 输出目录
-p 代表双端测序
--bam 代表输入文件为bam文件
-o 输出目录
3-2.提取甲基化信息（可选）
bismark_methylation_extractor -p --comprehensive --no_overlap --bedGraph --counts --buffer_size 20G --report --cytosine_report --genome_folder 基因组所在目录 -o 输出目录 输入的bam文件
-s 输入文件是从单端读取数据生成的Bismark结果文件。如果既没有-s也不-p设置实验的类型将被自动确定
-p 输入文件是从双端读取数据生成的Bismark结果文件。如果既没有-s也不-p设置实验的类型将被自动确定
--comprehensive 将可能的甲基化信息都输出，包括CHG,CHH,CpG
--no_overlap 对于双端读取，read_1和read_2理论上是可能重叠的。这个选项避免了重复的甲基化计算了2次。虽然这消除了对序列片段中心更多甲基化的计算偏差，但实际上可能会删除相当一部分数据，默认调用此参数
--bedGraph 输出bedGraph文件
--counts 指在bedGraph中有每个C上甲基化reads和非甲基化reads的数目(在--cytosine_report指定的情况下)
--buffer_size 使用指定的内存去计算甲基化信息
--report 会有一个甲基化的summary
--cytosine_report 指报道全基因组所有的CpG
--genome_folder 参考基因组的位置
-o 输出目录
重要文件有bismark_bt2_pe.CpG_report.txt（可以用来变化格式使methylkit包可以读取）和bismark_bt2_pe.bedGraph.gz（可以用IGV打开，报告给定胞嘧啶的位置及其甲基化状态）
3-3.总结甲基化信息（可选）
bismark2report（使用bismark_bt2_PE_report.txt）
bismark2summary（使用bismark_bt2_pe.bam）
4.排序
samtools sort 输入的bam文件 -o 输出的排序后的bam文件
(如要指定目录的话，目录一定要先创建好，其不能在运行时创建)
5.读入R
my.methRaw=processBismarkAln(location =bam文件全路径(用list表示),
sample.id=样品分组(用list表示),
assembly="mm10",# 比对时所用的参考基因组
read.context = "none",#要读入内存的内容
save.context=c("CpG","CHG","CHH"),#要保存的内容
save.folder=getwd())#保存路径
6.回到terminal，进入工作目录，去除非一般染色体
cat CpG文件 | grep chr > 新的CpG文件
7.重新读入R
myobj=methRead(location=CpG文件全路径(用list表示),
sample.id=样品分组(用list表示),
assembly="mm10",
treatment=c(1,1,0,0),
context="CpG",
mincov = 10)
示例：
library(methylKit)
file.list=list(system.file("extdata", "test1.myCpG.txt", package = "methylKit"),
system.file("extdata", "test2.myCpG.txt", package = "methylKit"),
system.file("extdata", "control1.myCpG.txt", package = "methylKit"),
system.file("extdata", "control2.myCpG.txt", package = "methylKit"))
myobj=methRead(file.list,
sample.id=list("test1","test2","ctrl1","ctrl2"),
assembly="hg38",
treatment=c(1,1,0,0),
context="CpG",
mincov = 10)
#理想情况下，您应该首先过滤具有极高（低）覆盖率的碱基，可使用filterByCoverage()函数解决PCR偏倚性，
#然后运行normalizeCoverage()函数以标准化样本之间的覆盖率。
#这两个功能将有助于减少统计测试中由于某些样本读数的系统过采样而可能产生的偏差。
getMethylationStats(myobj[[1]],plot=TRUE,both.strands=FALSE)#样品统计信息
getCoverageStats(myobj[[1]],plot=TRUE,both.strands=FALSE)#样品统计信息
myobj=filterByCoverage(myobj,lo.count=10,lo.perc=NULL,hi.count=NULL,hi.perc=99.9)
myobj <- normalizeCoverage(myobj)
tiles=tileMethylCounts(myobj,win.size=500,step.size=500)
meth=unite(tiles, destrand=FALSE)#合并样本_DMR
# meth=unite(myobj, destrand=FALSE)#合并样本_DML
head(meth)
#样品聚类信息
getCorrelation(meth,plot=TRUE)
clusterSamples(meth, dist="correlation", method="ward", plot=TRUE)
# hc = clusterSamples(meth, dist="correlation", method="ward", plot=F)
PCASamples(meth, screeplot=TRUE)
PCASamples(meth)
##查找差异甲基化碱基或区域
myDiff=calculateDiffMeth(meth)
myDiff25p=getMethylDiff(myDiff,difference=25,qvalue=0.01)
diffMethPerChr(myDiff,plot=T,qvalue.cutoff=0.01, meth.cutoff=25)
##统计甲基化位点的基因组分布和cpg岛分布(即注释差异甲基化的碱基或区域)
library(genomation)
myDiff_GRanges <- as(myDiff25p,"GRanges")
gene.obj=readTranscriptFeatures("c:/Users/ZSY/Desktop/DNA_methylation/hg38.bed.txt")#读取bed文件，有关启动子/外显子/内含子的注释
diffAnn=annotateWithGeneParts(myDiff_GRanges,gene.obj)# 差异甲基化位点的基因组分布
#library(TxDb.Hsapiens.UCSC.hg38.knownGene)
#diffAnn <- annotateDMRInfo(myDiff_GRanges, 'TxDb.Hsapiens.UCSC.hg38.knownGene')
cpg.obj=readFeatureFlank("c:/Users/ZSY/Desktop/DNA_methylation/hg38.bed.txt",feature.flank.name=c("CpGi","shores"))#读取CpG岛（CpGi）和CpG海岸注释
diffCpGann=annotateWithFeatureFlank(myDiff_GRanges,cpg.obj$CpGi,cpg.obj$shores,
feature.name="CpGi",
flank.name="shores")#差异甲基化位点的CpG岛分布
head(getMembers(diffAnn))#获得感兴趣的区域是否与外显子/内含子/启动子重叠
head(getMembers(diffCpGann))#获得感兴趣的区域是否与CpG岛重叠
head(getAssociationWithTSS(diffAnn))#获得距TSS的距离最近的基因名称
getTargetAnnotationStats(diffAnn,percentage=TRUE,precedence=TRUE)#与内含子/外显子/启动子重叠的差异甲基化区域的百分比/数目
plotTargetAnnotation(diffAnn,precedence=TRUE,main="differential methylation annotation")#上面的数据画图
getFeatsWithTargetsStats(diffAnn,percentage=TRUE)#与差异甲基化碱基重叠的内含子/外显子/启动子占所有的内含子/外显子/启动子的百分比
plotTargetAnnotation(diffCpGann,col=c("green","gray","white"),main="differential methylation annotation")#显示CpG岛，CpG海岸和其他地区差异甲基化碱基的百分比
上面所用的gene和cpg bed文件可从ucsc下载
