甲基化芯片数据的差异分析
测试开头
测试结尾
今天是生信星球陪你的第612天
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
花花写于2020-5-1，去年今天签了劳动合同，正式工作一周年啦~纪念一下。感谢曾老板！火山图有误，阅读原文看修改版
前情提要
前面两个步骤已经完成了数据下载、探针过滤、数据质控、样本过滤。
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)
插个小广告！
全国巡讲全球听（生信线上直播课）
请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！
