既然可以看感兴趣基因的生存情况，当然就可以批量做完全部基因的生存分析 / 开普饭

前言

就是一篇文章并没有使用TCGA数据库的指定癌症的生存信息去看自己感兴趣的基因的生存效应，反而舍近求远去下载BMC Cancer. 2011 文章数据，所以我怀疑TCGA应该是该基因在该癌症里面的生存效果不显著！

所以就安排学徒来完成，下面是他的表演：

接下来，对GSE20685的所有基因做生存分析（表达量中位值分组），获取统计学显著差异的基因列表。

1. 批量生存分析，获取统计学显著差异的基因列表

rm(list=ls())

options(stringsAsFactors = F)

options(warn = -1)
library(AnnoProbe)

gset<-geoChina("GSE20685")

eSet <- exprs(gset[[1]])

pheno <- pData(gset[[1]])

dim(eSet)

#[1] 54627   327

dim(pheno)

#[1] 327  63

eSet[1:4,1:4]

#          GSM519117 GSM519118 GSM519119 GSM519120

#1007_s_at 11.995510 12.042242 11.903921 11.905713

#1053_at    9.440330  9.329414  8.893978  9.577525

#117_at     7.907728  7.991689  8.163261  9.160542

#121_at     9.007045  8.976325  8.728571  8.665711

##表达矩阵行名是探针名，需要转换为基因名。

checkGPL(gset[[1]]@annotation) #[1] TRUE e2g <- idmap(gset[[1]]@annotation) eSet <- filterEM(eSet,e2g) eSet[1:4,1:4] # GSM519117 GSM519118 GSM519119 GSM519120 #ZZZ3 10.562788 10.68748 10.209871 10.48307 #ZZEF1 9.856933 10.46540 9.873843 10.07197 #ZYX 10.672518 10.48420 11.005133 10.59319 #ZYG11B 11.059366 11.29188 11.039351 11.49405

查看临床信息

characteristics_ch1.3 是存活情况
characteristics_ch1.4 是存活时间（存疑）

dat <- pheno[,c("characteristics_ch1.3","characteristics_ch1.4")] library(stringr) colnames(dat) <- c("status","OS.time") dat$status <- as.numeric(unlist(lapply(str_split(dat$status,":"),function(x) {return(x[2])}))) dat$OS.time <- as.numeric(unlist(lapply(str_split(dat$OS.time,":"),function(x) {return(x[2])}))) table(dat$status) # 0 1 #244 83 boxplot(dat$OS.time~dat$status)

死亡病例对应的持续回访时间更低。

library(survival)


cg <- array(dim = nrow(eSet))
survdata <- dat

my.surv <- Surv(survdata$OS.time,survdata$status==0)

#library("survminer")

dim(eSet)

for (i in 1:nrow(eSet)) { # i <- 1 # print(i) gene_exprs <- eSet[i,] gene_exprs <- as.data.frame(t(gene_exprs)) gene_exprs[,1] <- as.numeric(gene_exprs[,1]) gene_exprs$g <- ifelse(gene_exprs[,1]>=median(gene_exprs[,1]),'high','low') # print(table(gene_exprs$g)) gene_exprs$sample_name <- rownames(gene_exprs) gene_surv <- survdata gene_surv$sample_name <- rownames(gene_surv) gene_surv <- merge(gene_surv,gene_exprs,by='sample_name') # kmfit <- survfit(my.surv~gene_surv$g,data = gene_surv) kmdif <- survdiff(my.surv~gene_surv$g,data = gene_surv) # ggsurvplot(kmfit,conf.int = F,pval = T) p.value <- 1 - pchisq(kmdif$chisq, length(kmdif$n) -1) cg[i] <- (p.value < 0.05) } table(cg) #cg #FALSE TRUE #18386 1802 ##有1802个基因统计学差异显著。 surv_diff_genes <- rownames(eSet)[cg] save(surv_diff_genes,file = "surv_diff_genes.rdata")

2. 对生存分析显著的基因列表做富集分析

参考：为R包写一本书（像Y叔致敬）

获取列表基因的ENTREZID

rm(list=ls()) load("surv_diff_genes.rdata") surv.diff.genes <- as.data.frame(surv_diff_genes) colnames(surv.diff.genes) <- c("SYMBOL") library(clusterProfiler) library(org.Hs.eg.db) df <- bitr(surv.diff.genes$SYMBOL, fromType = "SYMBOL", toType = c("ENTREZID"), OrgDb = org.Hs.eg.db) genelist <- df$ENTREZID length(genelist) #[1] 1762

GO分析

go_enrich_results <- lapply( c('BP','MF','CC') , function(ont) { cat(paste('Now process ',ont )) ego <- enrichGO(gene = genelist, OrgDb = org.Hs.eg.db, ont = ont , pAdjustMethod = "BH", pvalueCutoff = 0.99, qvalueCutoff = 0.99, readable = TRUE)

print( head(ego) ) dotplot(ego,title=paste0('dotplot_',ont)) %>% print() return(ego) }) save(go_enrich_results,file = 'go_enrich_results.Rdata')

KEGG分析

library(clusterProfiler) kk <- enrichKEGG(gene = genelist, organism = 'hsa', pvalueCutoff = 0.9, qvalueCutoff = 0.9) head(kk)[,1:6] # ID Description GeneRatio BgRatio pvalue p.adjust #hsa04145 hsa04145 Phagosome 25/658 152/7978 0.0006218799 0.1940265 #hsa04974 hsa04974 Protein digestion and absorption 16/658 95/7978 0.0044579625 0.6954421 #hsa04657 hsa04657 IL-17 signaling pathway 15/658 94/7978 0.0095703403 0.8056927 #hsa00650 hsa00650 Butanoate metabolism 6/658 28/7978 0.0241405319 0.8056927 #hsa05130 hsa05130 Pathogenic Escherichia coli infection 25/658 202/7978 0.0259336441 0.8056927 #hsa03010 hsa03010 Ribosome 20/658 153/7978 0.0260288000 0.8056927 enrichKK=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')

分析结果可视化，经典五图流：

dotplot(enrichKK,color = "pvalue") barplot(enrichKK,showCategory=20,color = "pvalue") cnetplot(enrichKK, categorySize="pvalue", colorEdge = TRUE) emapplot(enrichKK,color = "pvalue") heatplot(enrichKK,showCategory = 20)

气泡图：