甲基化芯片数据的差异分析

火山图有误，已更新

前情提要

前面两个步骤已经完成了数据下载、探针过滤、数据质控、样本过滤。

rm(list=ls())
load("./Rdata/step2_filtered_pd_myNorm.Rdata")
dim(myNorm)
#> [1] 412481     52
dim(pd)
#> [1] 52  4

此处需要补充两个知识点：

beta值的生物学意义

beta>=0.6 完全甲基化
beta<=0.2 完全未甲基化
0.2<beta<0.6 部分甲基化

差异分析的三个等级

DMP 位点差异甲基化位点分析-limma
DMR差异甲基化区域分析（连续的差异片段）
DMB 更大的区域 /区域分类（某个基因附近的全部甲基化探针）

在这个例子里只做到DMP。

3. 差异分析

champ包非常强大，差异甲基化位点分析只用一个函数：champ.DMP完成。并且分析得到的结果数据里自带了注释，可以拿去做富集分析。


library(ChAMP)
library(tibble)
# 差异分析
group_list <- pd$group_list
myDMP <- champ.DMP(beta = myNorm,pheno=group_list)
#>          Contrasts
#> Levels    pTumor-pNormal
#>   pNormal             -1
#>   pTumor               1
head(myDMP$Tumor_to_Normal)
#>                logFC   AveExpr        t      P.Value    adj.P.Val        B
#> cg12825070 0.6503432 0.3999306 36.63169 3.303066e-38 1.362452e-32 76.65472
#> cg13912117 0.6118236 0.3402117 33.39330 3.032185e-36 6.253593e-31 72.29361
#> cg14416371 0.6663104 0.3610405 32.55210 1.047587e-35 1.421679e-30 71.09053
#> cg07176264 0.5591093 0.3749723 32.36840 1.378661e-35 1.421679e-30 70.82367
#> cg23690166 0.6019141 0.3553090 30.55516 2.242895e-34 1.850303e-29 68.10619
#> cg08089301 0.5219445 0.3114262 29.61805 1.005654e-33 6.509974e-29 66.63928
#>            Tumor_AVG Normal_AVG  deltaBeta CHR   MAPINFO Strand Type     gene
#> cg12825070 0.7251022 0.07475903 -0.6503432   5 148033708      F    I     HTR4
#> cg13912117 0.6461235 0.03429987 -0.6118236   8 132054555      R    I         
#> cg14416371 0.6941957 0.02788529 -0.6663104  11  43602847      R    I MIR129-2
#> cg07176264 0.6545270 0.09541771 -0.5591093   2 120281999      F    I     SCTR
#> cg23690166 0.6562660 0.05435192 -0.6019141  17  46711017      R    I MIR196A1
#> cg08089301 0.5723984 0.05045394 -0.5219445  17  46655561      F    I    HOXB4
#>            feature    cgi       feat.cgi    UCSC_CpG_Islands_Name DHS Enhancer
#> cg12825070 1stExon island 1stExon-island chr5:148033472-148034080  NA       NA
#> cg13912117     IGR island     IGR-island chr8:132052203-132054749  NA       NA
#> cg14416371  TSS200 island  TSS200-island  chr11:43602545-43603215  NA     TRUE
#> cg07176264 1stExon island 1stExon-island chr2:120281661-120282211  NA       NA
#> cg23690166 TSS1500 island TSS1500-island  chr17:46710812-46711419  NA       NA
#> cg08089301 1stExon island 1stExon-island  chr17:46655215-46655604  NA       NA
#>                                 Phantom Probe_SNPs Probe_SNPs_10
#> cg12825070                                                      
#> cg13912117                                                      
#> cg14416371                                                      
#> cg07176264 high-CpG:119998435-119998530  rs2244214              
#> cg23690166   high-CpG:44065942-44066037                         
#> cg08089301

df_DMP <- myDMP$Tumor_to_Normal
df_DMP=df_DMP[df_DMP$gene!="",]
logFC_t <- 0.45
P.Value_t <- 10^-15
df_DMP$change <- ifelse(df_DMP$adj.P.Val < P.Value_t & abs(df_DMP$logFC) > logFC_t,
                        ifelse(df_DMP$logFC > logFC_t ,'UP','DOWN'),'NOT') 
table(df_DMP$change) 
#> 
#>   DOWN    NOT     UP 
#>    345 108379    814
save(df_DMP,file = "./Rdata/step3.df_DMP.Rdata")

此处设置的logFC和p值的阈值是与原文一致的，由于甲基化的beta值不同于表达量，实际上用logFC也不是很对。曾老师课上推荐用deltabeta值替代logFC，就是甲基化信号值的差值。

火山图

library(dplyr)
library(ggplot2)
dat  = rownames_to_column(df_DMP)
for_label <- dat%>% head(3)
p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(adj.P.Val))) +
  geom_point(alpha=0.4, size=3.5, 
             aes(color=change)) +
  ylab("-log10(Pvalue)")+
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p

volcano_plot <- p +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(
    aes(label = rowname),
    data = for_label,
    color="black"
  )
volcano_plot

差异基因热图

cg <-  rownames(df_DMP[df_DMP$change != "NOT",])
plot_matrix <- myNorm[cg,]
annotation_col <- data.frame(Sample=pd$group_list) 
rownames(annotation_col) <- colnames(plot_matrix)
ann_colors = list(Sample = c(Normal="#4DAF4A", Tumor="#E41A1C"))

library(pheatmap)
pheatmap(plot_matrix,show_colnames = T,
         annotation_col = annotation_col,
         border_color=NA,
         color = colorRampPalette(colors = c("white","navy"))(50),
         annotation_colors = ann_colors,show_rownames = F)