代碼分享：基于TCGA甲基化數(shù)據(jù)挖掘驗(yàn)證分子表型代碼分享(第二期)

# a = column to rownames(methy sort,"V1")
# write.csv(me_data,'./00.data//pangguangai_me_data_412481.csv')
# write.csv(merge_info,'./00.data//merge_info.csv')
 
rm(list = ls())
options(stringsAsFactors = F)
library(ChAMP)
library(dplyr)
library(tibble)
library(openxlsx)
library(stringr)
library(data.table)
library(impute)
library(ggplot2)
library(ggpubr)
library(RColorBrewer)
library(pheatmap)
library(ggrepel)
set.seed(123)
 
getwd()
# setwd('C:/Users/30506/Desktop/CJ/易基因/研發(fā)部/CJ/公眾號(hào)文章復(fù)現(xiàn)/25.8.8-公共數(shù)據(jù)文章代碼分享-實(shí)操/數(shù)據(jù)處理文件夾/ DNA甲基化的探針和CG島的處理')
 
me_data <- read.csv("./pangguangai_me_data_412481.csv", header = T, row.names = 1)
merge_info<- read.csv("./merge_info.csv", header = T, row.names = 1)
 
beta_value = as.matrix(me_data)
#beta信號(hào)值矩陣?yán)锩娌荒苡蠳A值
# beta=impute.knn(beta_value)
# sum(is.na(beta))#這里是檢查矩陣?yán)锸欠襁�有NA
# beta=beta$data
# beta=beta+0.00001 #保證非0
#準(zhǔn)備pd表型文件(實(shí)際上就是樣品的信息)
pd_1<- as.data.frame(colnames(beta_value))
pd_info <- merge_info[merge_info$sample_submitter_id %in% pd_1$`colnames(beta_value)`,]
colnames(pd_1)="sample_submitter_id"
pd<- merge(pd_1,pd_info,by="sample_submitter_id",all.x= TRUE)#
colnames(beta_value)
colnames(merge_info)
 
myLoad=champ.filter(beta=beta_value,pd=pd)#beta就是指的是beta值，需要的甲基化信號(hào)
saveRDS(myLoad,file = './/pangguangai_21paired_meLoad_data_412481.rds')
 

# 5.質(zhì)控 -------------------------------------------------------------------
dir.create("../01.qc")
setwd('../01.qc')
QC =champ.QC(beta = myLoad$beta, pheno = myLoad$pd$sample_type)
#可以粗略的看出有些樣本需要去掉
#
myNorm<- champ.norm(beta=myLoad$beta,arraytype="450K",cores=4)
dim(myNorm)
 
#歸一化后會(huì)產(chǎn)生很多NA,看一下哪些樣品里含有NA,并去掉
num.na<- apply(myNorm,2,function(x)(sum(is.na(x))))
table(num.na)
names(num.na)= colnames(myNorm)
dt = names(num.na[num.na>0])
dn = str_replace(dt,"-01","-11")
keep=setdiff(colnames(myNorm),c(dt,dn))#只取colnames(myNorm)里的元素
myNorm = myNorm[,keep]
dim(myNorm)
pd = myLoad$pd
pd <- pd[pd$sample_submitter_id %in% colnames(myNorm),]
dim(pd)
 
##去掉對(duì)應(yīng)的離群樣本，下面的三個(gè)樣本的PCA不在一起
exp_count_gene_name =myLoad$beta
sample= pd
normalized_counts =myNorm
 
exp_count_gene_name = exp_count_gene_name[,c(1:10,13:18,21:24,27:30,33:42)]
sample=sample[c(1:10,13:18,21:24,27:30,33:42),]
normalized_counts=normalized_counts[,c(1:10,13:18,21:24,27:30,33:42)]
 
color_up <- colorRampPalette(c(rev(colorRampPalette(brewer.pal(11, "RdBu"))(1100))[701:1100]))(100)
color_middle <- colorRampPalette(c(rev(colorRampPalette(brewer.pal(11, "RdBu"))(1100))[401:700]))(100)
color_down <- colorRampPalette(c(rev(colorRampPalette(brewer.pal(11, "RdBu"))(1100))[1:400]))(100)
color_bar <- c(color_down, color_middle, color_up)
# colorRampPalette(brewer.pal(8, "Reds"))(100)
 
# setwd('./04.Comb.Count.PCA.250829/')
# normalized_counts =myNorm
# colnames(normalized_counts)
# head(myLoad$beta)
# exp_count_gene_name =myLoad$beta
# colnames(exp_count_gene_name)
 
cor.pearson <- cor(exp_count_gene_name, exp_count_gene_name, method = "pearson")
dir.create('../01.Cor_Heat')
pdf("../01.Cor_Heat/01.1.Cor_Count.me.CJ.2508.pdf", width = 25, height = 25)
pheatmap(cor.pearson, border_color = NA,
         color = color_bar,
         cluster_cols = T, cluster_rows = T,
         display_numbers = TRUE, fontsize_number = 8, number_format = "%.4f",
         # fontsize=9, fontsize_row=6,
         # scale = 'row',
         show_rownames = T, show_colnames = T,
         legend = T,
         cellwidth = 30, cellheight = 30)
dev.off()
 
cor.pearson <- cor(log2(exp_count_gene_name+1), log2(exp_count_gene_name+1), method = "pearson")
 
pdf("../01.Cor_Heat/01.1.Cor_log2Count.me.CJ.2508.pdf", width = 25, height = 25)
pheatmap(cor.pearson, border_color = NA,
         color = color_bar,
         cluster_cols = T, cluster_rows = T,
         display_numbers = TRUE, fontsize_number = 8, number_format = "%.4f",
         # fontsize=9, fontsize_row=6,
         # scale = 'row',
         show_rownames = T, show_colnames = T,
         legend = T,
         cellwidth = 30, cellheight = 30)
dev.off()
 
# dir.create('../02.PCA')
project.pca <- prcomp(t(log2(normalized_counts+1)))
# project.pca <- prcomp(t(normalized_counts))
summary(project.pca)
 
#Determine the proportion of variance of each component
#Proportion of variance equals (PC stdev^2) / (sum all PCs stdev^2)
project.pca.proportionvariances <- ((project.pca$sdev^2) / (sum(project.pca$sdev^2)))*100
 
barplot(project.pca.proportionvariances, cex.names=1, xlab=paste("Principal component (PC), 1-", length(project.pca$sdev)), ylab="Proportion of variation (%)", main="Scree plot", ylim=c(0,100))
# ggsave("../02.PCA/barplot_PCA.pdf", width = 7, height = 6)

# ************************ caculate by prcomp ************************ #
# sample= pd
pcaData.prcomp <- data.frame(row.names = colnames(normalized_counts),
                             PC1 = project.pca$x[,1],
                             PC2 = project.pca$x[,2],
                             PC3 = project.pca$x[,3],
                             PC4 = project.pca$x[,4],
                             PC5 = project.pca$x[,5],
                             Condition = sample$case_submitter_id,
                             Replicate = sample$sample_type)
 
color <- colorRampPalette(brewer.pal(8, "Set1"))(13) # !!!挑的這個(gè)
# 互補(bǔ)色 30° * 12
color <- c("#FF7373", "#FF0000",
           "#FFB273", "#FF7400",
           "#FFD073", "#FFAA00",
           "#FFE773", "#FFD300",
           "#FFFF73", "#FFFF00",
           "#C9F76F", "#9FEE00",
           "#67E667", "#00CC00",
           "#5CCCCC", "#009999",
           "#6C8CD5", "#1240AB",
           "#876ED7", "#3914AF",
           "#AD66D5", "#7109AA",
           "#E667AF", "#CD0074"
)
 
# color <- c("#FF7373", "#FF0000",
#            "#FFC373", "#FF9200",
#            "#FFE773", "#FFD300",
#            "#E3FB71", "#CCF600",
#            "#67E667", "#00CC00",
#            "#66A3D2", "#0B61A4",
#            "#876ED7", "#3914AF",
#            "#D25FD2", "#A600A6"
# )
 
p1 <- ggplot(pcaData.prcomp, aes(PC1, PC2, color=Condition, shape = Replicate)) +
  geom_point(size=4.5) +
  scale_color_manual(values = color) +                           # scale_color_manual 手動(dòng)設(shè)置顏色, color對(duì)color，fill對(duì)fill
  xlab(paste0("PC1: ", round(project.pca.proportionvariances[1], 0), "% variance")) +
  ylab(paste0("PC2: ", round(project.pca.proportionvariances[2], 0), "% variance")) +
  # xlim(-200,200) +
  # ylim(-200,200) +
  coord_fixed() +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(color = "black", fill = NA))                   # 去除背景和網(wǎng)格線，并加上黑色邊框
p1
# dir.create('./02.PCA')   #推薦PC1  PC2
ggsave("../02.PCA/01.PC1_PC2_log2Count.Batch2508.pdf", width = 7, height = 6)
# ggsave("02.PCA/01.PC1_PC3_log2Count.Batch2508.pdf", width = 7, height = 6)
# ggsave("02.PCA/01.PC2_PC3_log2Count.Batch2508.pdf", width = 7, height = 6)

 
p2 <- p1 + geom_text_repel(aes(x = PC1, y = PC2,
                               label = colnames(normalized_counts)),
                           colour="darkred", size = 2, box.padding = unit(0.35, "lines"),
                           point.padding = unit(0.3, "lines"))
p2
ggsave("../02.PCA/01.PC1_PC2_log2Count.Label.Batch2508.pdf", width = 14, height = 12)

dir.create('../03.plot')
n=names(tail(sort(apply(normalized_counts,1,sd)),500))
 
ac=data.frame(group=sample$sample_type)
rownames(ac)=colnames(normalized_counts)
color_bar <- colorRampPalette(brewer.pal(9, "Purples"))(100)
pheatmap(normalized_counts[n,],show_colnames =F,show_rownames = F, cluster_row = T,
         annotation_col=ac,treeheight_row = 0,treeheight_col = 0,color = color_bar,
         filename = '../03.plot/heatmap_normalized_counts.pdf'
)
 
# dev.off()
color_bar <- colorRampPalette(brewer.pal(9, "Purples"))(100)
color_bar <-colorRampPalette(c("#89B9D7","white","#C92A1F"))(50)
pheatmap(normalized_counts[n,],show_colnames =F,show_rownames = F, cluster_row = T,
         annotation_col=ac,treeheight_row = 0,treeheight_col = 0,color = color_bar,
         filename = '../03.plot/heatmap_normalized_counts_1.pdf'
)

color_bar <-colorRampPalette(c("#625D9E","white","#B53E2B"))(50)
pheatmap(normalized_counts[n,],show_colnames =F,show_rownames = F, cluster_row = T,
         annotation_col=ac,treeheight_row = 0,treeheight_col = 0,color = color_bar,
         filename = '../03.plot/heatmap_normalized_counts_2.pdf'
)
 

# color_bar <-colorRampPalette(c("#625D9E","white","#B53E2B"))(50)
pheatmap(normalized_counts[n,],show_colnames =F,show_rownames = F, cluster_row = T,
         annotation_col=ac,treeheight_row = 0,treeheight_col = 0,
         filename = '../03.plot/heatmap_normalized_counts_3.pdf'
)
#上面是挑選的500個(gè)sd最大的cg

# 6.差異分析 ------------------------------------------------------------------
#DMPs 分析
sample$sample_type_fix = sample$sample_type
# sample$sample_type_fix
sample$sample_type_fix[sample$sample_type_fix=='Solid Tissue Normal']='Normal'
sample$sample_type_fix[sample$sample_type_fix=='Primary Tumor']='Tumor'
 
myDMP <- champ.DMP(beta = normalized_counts,pheno=sample$sample_type_fix)
#查看結(jié)果
head(myDMP[[1]])
 
df_DMP <- myDMP$Tumor_to_Normal
df_DMP=df_DMP[df_DMP$gene!="",]
head(df_DMP)
#參數(shù)需要調(diào)整為合適
logFC_t <- 0.4
P.Value_t <- 10^-5
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'),'None')
table(df_DMP$change)
# Down  None    Up
# 1195 80244  1115
dat =rownames_to_column(df_DMP)
for_label <- dat%>% head(3)
mycol <- c("#EB4232","#2DB2EB","#d8d8d8")
p1=ggplot(data = dat,aes(x = logFC,y=-log10(adj.P.Val)))+
  geom_point(alpha=0.4,size=3.5,aes(color=change))+
  ylab("-log10(adj.P.Val)")+
  scale_color_manual(values=c("green","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()
p1
ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5.pdf"), width = 6*1.2, height = 6*1)
ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5.png"), width = 6*1.2, height = 6*1)
 

# dev.off()
# dat$change
p <- ggplot(data = dat,
            aes(x = logFC, y = -log10(adj.P.Val), color = change)) + #建立映射
  geom_point(size = 2.2) #繪制散點(diǎn)
p
 
#自定義顏色與主題：
mycol <- c("#2DB2EB","#d8d8d8","#EB4232")
 
mytheme <- theme_classic() +
  theme(
    axis.title = element_text(size = 15),
    axis.text = element_text(size = 14),
    legend.text = element_text(size = 14)
  )
p2 <- p +
  scale_colour_manual(name = "", values = alpha(mycol, 0.7)) +
  mytheme
p2
 
#添加閾值線：
p3 <- p2 +
  geom_hline(yintercept = -log10(P.Value_t),size = 0.7,color = "black",lty = "dashed") + #水平閾值線
  geom_vline(xintercept = c(-logFC_t,logFC_t),size = 0.7,color = "black",lty = "dashed") #垂直閾值線
p3
 
##在上調(diào)和下調(diào)里面找到對(duì)應(yīng)的cg
marker = read.xlsx('../00.data/副本Promoter_methylation_candidate(1).xlsx')
marker=marker[1:10,]
marker = dat[dat$gene %in% marker$SYMBOL,]
marker = marker[marker$change %in% c('Up','Down'),]
# marker$gene
p7 <- p3 +
  geom_text_repel(data = marker,
                  aes(x = logFC, y = -log10(adj.P.Val), label = gene),
                  force = 80, color = '#B53E2B', size = 5,
                  point.padding = 0.5, hjust = 0.1,
                  # arrow = arrow(length = unit(0.02, "npc"),
                  #               type = "open", ends = "last"),
                  segment.color="black",
                  segment.size = 0.3,
                  nudge_x = 0,
                  nudge_y = 1
  )
p7
ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5_1.pdf"), width = 6*1.2, height = 6*1)
ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5_1.png"), width = 6*1.2, height = 6*1)
 

marker$rowname
p8 <- p3 +
  geom_text_repel(data = marker,
                  aes(x = logFC, y = -log10(adj.P.Val), label = rowname),
                  force = 80, color = '#B53E2B', size = 5,
                  point.padding = 0.5, hjust = 0.1,
                  # arrow = arrow(length = unit(0.02, "npc"),
                  #               type = "open", ends = "last"),
                  segment.color="black",
                  segment.size = 0.3,
                  nudge_x = 0,
                  nudge_y = 1
  )
p8
ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5_1_cg.pdf"), width = 6*1.2, height = 6*1)
ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5_1_cg.png"), width = 6*1.2, height = 6*1)
 

save(merge_info,myDMP,myNorm,marker,df_DMP,
     file = "../00.data//pangguangai_21paired_meLoad_data_412481.Rdata")
 
#甲基化的信號(hào)值可以是β—value
head(hm450.manifest.hg19)
load('C:/Users/chang/Desktop/CJ /hm450_hg38_manifest_new-from-findOL.rda')
head(myDMP)
##把hg19的chrom換成這里的對(duì)應(yīng)pos的位置
cd=intersect(man38new$cgid,rownames(df_DMP))
dim(df_DMP) #82554
length(cd) #80573  少了2000多個(gè)探針位點(diǎn)對(duì)應(yīng)不到
pos=match(rownames(df_DMP),man38new$cgid,)
man38new_me=man38new[pos,]
man38new_me#
# man38new_me <- na.omit(man38new_me) #80573
colnames(df_DMP)
df_DMP$MAPINFO_hg38 = man38new_me$pos
 
df_DMP$MAPINFO_hg38
save(merge_info,myDMP,myNorm,marker,df_DMP,
     file = "../00.data//pangguangai_21paired_meLoad_data_412481.Rdata")
head(man38new) ##根據(jù)這個(gè)轉(zhuǎn)換過(guò)的文件對(duì)之前的代碼區(qū)域進(jìn)行替換即可
write.csv(df_DMP,'../00.data//df_DMP.csv')

# df_DMP
df_DMP_up =df_DMP[df_DMP$change %in% c('Up'),]
 
ego <- enrichGO(gene = unique(df_DMP_up$gene), OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="BP", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
 
go=data.frame(
  pathNames=ego@result$Description,
  count=ego@result$Count,
  Pvalue=ego@result$pvalue,
  P.adj=ego@result$p.adjust,
  log10_P_value <- -log10(ego@result$pvalue),##新增一列-log10pvalue
  geneID = ego@result$geneID,
  Hit.Ratio=ego@result$GeneRatio##該通路上你的基因的占比
)
colnames(go)[5] = 'log10_P_value'
dir.create('../03.plot/00.go')
write.csv(go,'../03.plot/00.go/up_gene_go.csv')
 
 
##篩選一下
ad =go[rownames(go) %in% c(854,1044,1626,1770,270,
                           841,1203,848,182,757
),]
go =go[go$pathNames %in% c(ad$pathNames),]
# library(paletteer)
Fig_1 <-
  go  %>%
  ggplot(aes(x = reorder(pathNames, log10_P_value), y = log10_P_value, label =  paste0("count=", signif(count)),
             fill ='up_gene_Bladder_cancer_Top10')) +
  theme_pubr(base_size = 20)+
  theme(axis.title.y = element_blank(),legend.title = element_blank(),
        axis.text.y = element_text(size = 15), axis.text.x = element_text(size = 18),
        legend.position = "top", strip.text = element_blank(), legend.key.size = unit(10, "pt"))+
  # labs(title = "WT_F-PGCs")+
  geom_bar(stat = "identity", alpha = 0.75)+
  # facet_grid(Cluster~., scales = "free", space = "free") +
  geom_text(hjust = 1.1, colour = "white", size = 3) +
  scale_fill_manual(values = '#E39A35') +
  coord_flip()
Fig_1
ggsave(paste0("../03.plot/00.go/", "up_gene_go_select",".pdf"),width = 6*2*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/00.go/", "up_gene_go_select",".png"),width = 6*2*1.5, height = 6*1.5)

df_DMP_dn =df_DMP[df_DMP$change %in% c('Down'),]
 
ego <- enrichGO(gene = unique(df_DMP_dn$gene), OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="BP", pvalueCutoff = 0.05, qvalueCutoff = 0.05)
 
go=data.frame(
  pathNames=ego@result$Description,
  count=ego@result$Count,
  Pvalue=ego@result$pvalue,
  P.adj=ego@result$p.adjust,
  log10_P_value <- -log10(ego@result$pvalue),##新增一列-log10pvalue
  geneID = ego@result$geneID,
  Hit.Ratio=ego@result$GeneRatio##該通路上你的基因的占比
)
colnames(go)[5] = 'log10_P_value'
dir.create('../03.plot/00.go')
write.csv(go,'../03.plot/00.go/dn_gene_go.csv')
 
##篩選一下
ad =go[rownames(go) %in% c(4:13),]
 
go =go[go$pathNames %in% c(ad$pathNames),]
# library(paletteer)
Fig_1 <-
  go  %>%
  ggplot(aes(x = reorder(pathNames, log10_P_value), y = log10_P_value, label =  paste0("count=", signif(count)),
             fill ='dn_gene_Bladder_cancer_Top10')) +
  theme_pubr(base_size = 20)+
  theme(axis.title.y = element_blank(),legend.title = element_blank(),
        axis.text.y = element_text(size = 15), axis.text.x = element_text(size = 18),
        legend.position = "top", strip.text = element_blank(), legend.key.size = unit(10, "pt"))+
  # labs(title = "WT_F-PGCs")+
  geom_bar(stat = "identity", alpha = 0.75)+
  # facet_grid(Cluster~., scales = "free", space = "free") +
  geom_text(hjust = 1.1, colour = "white", size = 3) +
  scale_fill_manual(values = '#57C3F3') +
  coord_flip()
Fig_1
ggsave(paste0("../03.plot/00.go/", "dn_gene_go_select",".pdf"),width = 6*2*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/00.go/", "dn_gene_go_select",".png"),width = 6*2*1.5, height = 6*1.5)