对于如何使用minfi 分析甲基化芯片数据，我们在之前的文章中详细讲解了每一步处理的具体用法。今天主要给出一个piepeline, 包括从文件读取一直到最终的DMP/DMR差异结果。
整个pipeline 可以分成以下几步：
- 数据读取
- 质量过滤
- 预处理
- 探针筛选
- 差异分析
完整的代码如下：
代码语言：javascript复制# 加载 minfi 包
library(minfi)# 根据SampleSheet.csv 文件，读取所有样本的 .idat 文件
targets <- read.metharray.sheet("./", pattern="HumanMethylation450_Demo_Sample_Sheet.csv")
rgSet <- read.metharray.exp(targets=targets)# 计算探针的p值，过滤掉在任何以一个样本中p值大于0.01的探针
probeP <- detectionP(rgSet)
keep <- apply(probeP, 1 , function(t){all(t < 0.01)})
rgSet <- rgSet[keep,]# 过滤掉所有探针p值的均值大于0.05的样本
keep <- apply(probeP, 2, mean) < 0.05
rgSet <- rgSet[,keep]# 预处理，背景降噪和归一化，注意，此处可以根据情况，替换成另外的算法
Gset.funnorm <- preprocessFunnorm(rgSet)# 探针过滤，去除在性染色体上的探针
annotation <- getAnnotation(Gset.funnorm)
sex_probe <- rownames(annotation)[annotation$chr %in% c("chrX", "chrY")]
keep <- !(featureNames(Gset.funnorm) %in% sex_probe)
Gset.funnorm <- Gset.funnorm[keep,]# 去除覆盖了SNP位点的探针，如果感觉过滤掉的探针太多，可以适当上调SNP频率, 将maf的值变大
GRset <- dropLociWithSnps(Gset.funnorm, snps=c("SBE","CpG"), maf=0)# DMP, 探针水平的差异分析
beta <- getBeta(GRset)
group <- pData(GRset)$Sample_Group
dmp <- dmpFinder(beta, pheno = group , type = "categorical")
head(dmp)# DMR， 甲基化区域的差异分析
group <- pData(GRset)$Sample_Group
designMatrix <- model.matrix(~ group)
dmrs <- bumphunter(GRset,
design = designMatrix,
cutoff = 0.2, B=0, type="Beta")
head(dmrs$table[,1:4], n =3 )总结
通过整个pipeline的代码，可以看出分析其实并不复杂，关键是要理解每个步骤处理的意义和结果的解读。这里的代码并没有给出一些可视化的图表结果，其实在minfi中也提供了可视化的解决方案，后面我们再详细探索。
