minfi 分析甲基化芯片数据-pipeline篇(附完整代码)

生信修炼手册

发布于 2020-05-12 10:50:18

1.8K0

文章被收录于专栏：生信修炼手册生信修炼手册

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

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

数据读取
质量过滤
预处理
探针筛选
差异分析

完整的代码如下：

# 加载 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 )