对于如何使用minfi 分析甲基化芯片数据，我们在之前的文章中详细讲解了每一步处理的具体用法。今天主要给出一个piepeline, 包括从文件读取一直到最终的DMP/DMR差异结果。

整个pipeline 可以分成以下几步：

1. 数据读取
2. 质量过滤
3. 预处理
4. 探针筛选
5. 差异分析

完整的代码如下：

# 加载 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中也提供了可视化的解决方案，后面我们再详细探索。