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]]