TCGA甲基化芯片数据的差异基因的生存分析与富集分析

4.将差异甲基化位点拿来做批量生存分析

rm(list = ls())
library(data.table)
library(stringr)
library(survival)
library(survminer)
load("./Rdata/step2_filtered_pd_myNorm.Rdata")
load("./Rdata/step3.df_DMP.Rdata")

cg <- rownames(df_DMP[df_DMP$change != "NOT",])
myNorm_tumor <- myNorm[cg,]

logrank test批量生存分析

这里和TCGA系列里的内容是一样的，详见：两种方法批量做生存分析

suv_dat <- data.table::fread("./raw_data/TCGA-HNSC.survival.tsv.gz")
suv_dat$sample = str_sub(suv_dat$sample,1,15)
suv_dat <- suv_dat[suv_dat$sample %in% colnames(myNorm_tumor),]
suv_dat <- suv_dat[str_sub(suv_dat$sample,14,15)=="01",] 
suv_dat = merge(pd,suv_dat,by.x = "sampleID",by.y = "sample")
myNorm_tumor <- myNorm_tumor[,suv_dat$sample]
identical(colnames(myNorm_tumor),suv_dat$sample)
#> [1] TRUE
library(survival)
logrankP <- apply(myNorm_tumor, 1, function(x){
  #x <- myNorm_tumor[1,]
  suv_dat$group <- ifelse(x>mean(x),"High","Low")
  res <- coxph(Surv(OS.time, OS)~group, data=suv_dat)
  beta <- coef(res)
  se <- sqrt(diag(vcov(res)))
  p.val <- 1 - pchisq((beta/se)^2, 1)
})
table(logrankP<0.05) #110个CpG位点
#> 
#> FALSE  TRUE 
#>  1676   110

得到110个对生存影响显著的差异甲基化位点，取前20个画图。

surv_gene <- names(sort(logrankP))[1:20] 
choose_matrix <- myNorm[surv_gene,]
annotation_col <- data.frame(Sample=pd$group_list) 
rownames(annotation_col) <- colnames(choose_matrix)
ann_colors = list(Sample = c(Normal="#4DAF4A", Tumor="#E41A1C"))

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

image.png

也可以画画生存分析的图

gs=head(surv_gene,4)
exprSet = myNorm_tumor
meta = suv_dat
splots <- lapply(gs, function(g){
  meta$gene=ifelse(exprSet[g,]>median(exprSet[g,]),'high','low')
  sfit1=survfit(Surv(OS.time, OS)~gene, data=meta)
  ggsurvplot(sfit1,pval =TRUE, data = meta)
}) 
arrange_ggsurvplots(splots, print = TRUE,  
                    ncol = 2, nrow = 2)

image.png

5.富集分析

利用ChAMP包对过滤后的数据做了差异甲基化位点分析。如果是肿瘤数据的话，可以加一步生存分析。

rm(list = ls())
load(file = 'Rdata/step3.df_DMP.Rdata')
library(ggplot2)
library(stringr)
library(clusterProfiler)
library(org.Hs.eg.db)

ID转换

length(unique(df_DMP$gene))
#> [1] 16338
s2e <- bitr(unique(df_DMP$gene), fromType = "SYMBOL",
           toType = c( "ENTREZID"),
           OrgDb = org.Hs.eg.db)

df_DMP=merge(df_DMP,s2e,by.y='SYMBOL',by.x='gene')
table(!duplicated(df_DMP$ENTREZID))
#> 
#> FALSE  TRUE 
#> 83778 14188

gene_up= unique(df_DMP[df_DMP$change == 'UP','ENTREZID'] )
gene_down=unique(df_DMP[df_DMP$change == 'DOWN','ENTREZID'] )
gene_diff=c(gene_up,gene_down)
gene_all=unique(df_DMP$ENTREZID)

富集分析及其可视化

之前完整的富集分析会把上下调基因分开、合并都做。GO富集分析会将MF、CC、BP三个部分分开做。代码太多了，在这简化一下。

kkgo_file = "./Rdata/kkgo_file.Rdata"
if(!file.exists(kkgo_file)){
  kk <- enrichKEGG(gene         = gene_diff,
                   universe     = gene_all,
                   organism     = 'hsa',
                   pvalueCutoff = 0.05)
  go <- enrichGO(gene_diff, OrgDb = "org.Hs.eg.db", ont="all") 
  save(kk,go,file = kkgo_file)
}
load(kkgo_file)

dotplot(kk)

image.png

barplot(go, split="ONTOLOGY",font.size =10)+ 
  facet_grid(ONTOLOGY~., scale="free") + 
  scale_x_discrete(labels=function(x) str_wrap(x, width=45))

image.png