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

 今天是生信星球陪你的第614天

---

   大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

   就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

   这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

接着差异分析往下做生存分析和富集分析~

当然，生存分析只是肿瘤数据可以做，需要对应的病人生存信息。不是肿瘤的数据就直接跳过生存分析  。目录等下次发咯~

大家有没有发现，蓝色排版的那个豆豆消失好几天了，请大家一起到文末投票，召唤豆豆。

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)

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))