methylKit：差异甲基化分析

1. 1Methylation输入格式

coverage：覆盖这个位点的reads数目。
freqC：甲基化C的比例
freqT：非甲基化C的比例

1.2 纯文本的读取：methRead

file.list是一个列表，里面的元素是每个文件名
myobj = methRead(file.list,
			   	 assembly = "hg19",
			   	 sample.id = list("test1","test2","ctrl1","ctrl2"),
			   	 treatment = c(1,1,0,0),
			   	 context = "CpG"
			   	 )
treatment指定样本的分组，1是treatment组，0是control组；assembly指定组装版本。

1.3 bam文件读取：processBismarkAln

bam.list是一个列表，里面元素是每个bam文件名。
myobj = processBismarkAln(bam.list,
						  sample.id = c("sample_A1","sample_A2","sample_B1","sample_B2"),
						  assembly = "hg19",
						  treatment = c(1,1,0,0),
						  read.context = "CpG",
						  save.fold = getwd())
						  

2.1 探索样本的描述性数据：getMethylationStats
myobj是一个methylRawList对象，它记录了每个样本的甲基化信息。

getMethylationStats(myobj[[2]],plot = FALSE/TRUE,both.strands = FALSE)
获取第二个样本的甲基化信息。

同时可以读取myobj测序深度信息。
getCoverageStats(myobj[[2]],plot = TRUE,both.strands = FALSE)

2.2 根据测序深度过滤样本：filterByCoverage

filtered.myobj = filterByCoverage(myobj,
							      lo.count = 10,
							      lo.perc = NULL,
							      hi.count = NULL,
							      hi.perc = 99.9
							      )
lo.count:指定最小深度
hi.count:指定最大深度
lo.perc:
hi.perc:

3 合并样本数据：unite
meth = unite(myobj,destrand = FALSE)

4.1样本的相关性：getCorrelation
getCorrelation(meth,plot = TRUE)

4.2样本聚类：clusterSamples、PCASamples
clusterSamples(meth,dist = "correlation",method = "ward",plot = TRUE)

PCA碎石图：
PCASamples(meth,screeplot= TRUE)

PCA主成分图：
PCASamples(meth)

5 批量效应
批量效应用于检查哪些主要成分与潜在的批处理效果上统计相关。如果某些主要成分是批次效应造成的结果，就需要使用removeComp将其删除。

sampleAnnotation = data.frame(batch_id = c("a","a","b","b"),
						      age = c(19,34,23,40))
as = assocComp(mBase = meth,sampleAnnotaion)
移除第一个主成分：
newObj=removeComp(meth,comp=1)

6 平铺窗口分析：tileMethylCounts
差异甲基化区域的分析:按照滑动窗口的方式定义甲基化区域，默认窗口大小为10000 bp ，步长为10000bp.
tiles=tileMethylCounts(myobj,win.size=1000,step.size=1000)

7 查找差异甲基化的碱基或区域
测试差异甲基化：

myDiff = calculateDiffMeth(meth)

基于q值、甲基化差异阈值百分比和type类型选择差异甲基化区域：

myDiff25p = getMethylDiff(myDiff,difference = 25,qvalue = 0.01,type = "hyper" / "hypo" / "all")

hyper表示相比control组，treatment组中的甲基化C更多；
hypo则相反，表示treatment组中的甲基化C比control组中少。

8 注释差异甲基化的碱基或区域：genomation包
读取基因BED文件：
gene.obj = readTranscriptFeatures("gene.bed")
使用BED文件注释甲基化区域：启动子、外显子、内含子
diffAnn = annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj)

注释后，获取距离TSS距离最近的基因名称：
getAssociationWithTSS(diffAnn)

区域分析例如汇总启动子区域集上的甲基化信息：

myobj：一开始读入的methylation信息。
gene.obj：读入的BED文件。
promoters = regionCounts(myobj,gene.obj$promoters)
promoters[[1]]