学习资料¶

DNA甲基化芯片分析视频：【生信技能树】DNA甲基化芯片分析
视频中完整代码见：https://github.com/jmzeng1314/methy_array
甲基化的一些基础知识：点击这里
450甲基化基础知识：https://www.jianshu.com/p/f7428021de9c

甲基化芯片分析流程¶

甲基化芯片数据介绍：

注：以下只展示了部分代码，完整代码见：https://github.com/jmzeng1314/methy_array

1、读取甲基化芯片数据（甲基化信号矩阵、idat文件）¶

1.1 读取甲基化信号矩阵 (也可以直接从GEO官网下载 GSE68777)¶

require(GEOquery)    
GSE68777 <- getGEO("GSE68777",destdir = './',AnnotGPL = T,getGPL = F) # 下载数据
# 获取表达矩阵
beta.m <- exprs(eset[[1]])
# 获取注释信息
pD.all <- pData(eset[[1]])

1.2 导入idat文件¶

使用minfi包

# 使用minfi
library("minfi")
rgSet <- read.metharray.exp("GSE68777/idat") # idat 文件存放目录
# 返回对象 class: RGChannelSet

使用ChAMP包

ChAMP 提供了完整的分析甲基化芯片的pipieline, 从数据导入到差异分析和下游的基因功能分析

# BiocManager::install("ChAMP",ask = F,update = F)
library("ChAMP")
# 第一个参数为原始的.idat文件和SampleSheet.csv 文件存放的目录，第二个参数为芯片的类型，包括450K和EPIC两种，默认为450K
myLoad <- champ.load(testDir,arraytype="450K")

2、质量检测¶

filter

library(ChAMP)
# beta 信号值矩阵里面不能有NA值
myLoad = champ.filter(beta = beta ,pd = pd)
# 标准化
myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
dim(myNorm) 
pD=myLoad$pd

dat=t(beta.mdat[1:4,1:4]) # dat.adat 是处理好的beta矩阵
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")  
# 因为甲基化芯片是450K或者850K，几十万行的甲基化位点，所以PCA不会太快
dat.pca <- PCA(dat , graph = FALSE) 
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = group_list, # color by groups
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)
ggsave('all_samples_PCA.png')

dat=beta.m
dat[1:4,1:4] 
cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=group_list)
rownames(ac)=colnames(n)  
pheatmap(n,show_colnames =F,show_rownames = F,
         annotation_col=ac,filename = 'heatmap_top1000_sd.png')
dev.off()

exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
dim(exprSet)
# M=cor(log2(exprSet+1)) 
M=cor(exprSet)
pheatmap::pheatmap(M,annotation_col = colD)
pheatmap::pheatmap(M,
                   show_rownames = F,
                   annotation_col = colD,
                   filename = 'cor_top500.png')
dev.off()

fig2

fig3

fig5

3、差异分析¶

ChAMP包

myDMP <- champ.DMP(beta = myNorm, pheno = group_list)

minfi包

dmp <- dmpFinder(M, pheno=group_list, type="categorical")
dmpDiff=dmp[(dmp$qval<0.05) & (is.na(dmp$qval)==F),]

4、GO或者KEGG等数据库的功能注释¶

load(file = 'step3-output-myDMP.Rdata')
deg=myDMP[[1]]
head(deg)
length(unique(deg$gene)) 
deg$g=ifelse(abs(deg$logFC) < 0.2,'stable',
         ifelse(deg$logFC > 0.2,'UP','DOWN'))

deg$symbol=deg$gene

library(org.Hs.eg.db)
df <- bitr(unique(deg$symbol), fromType = "SYMBOL",
           toType = c( "ENTREZID"),
           OrgDb = org.Hs.eg.db)
head(df)
DEG=deg
head(DEG)

DEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')