其他
既然可以看感兴趣基因的生存情况，当然就可以批量做完全部基因的生存分析
接下来，对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叔致敬）
01获取列表基因的ENTREZIDrm(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
02GO分析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')
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)
气泡图：
条带图：
通路与基因之间的关系可视化：
有趣的是生存分析有统计学显著的这些基因，富集得到的通路直接共享基因很少！
通路与通路之间的连接展示：
热图展现通路与基因之间的关系：
跑KEGG分析的时候，代码没变，跑了两次，出了两个结果。第一次只富集到1条通路，第二次跑的结果比较正常，2.2节所示。
估计第一次跑KEGG分析的时候，enrichKEGG函数下载数据的时候出现了问题。所以总结经验就是以后跑KEGG如果只跑出来1条通路，第一考虑enrichKEGG函数，重启R再跑一遍enrichKEGG。
生存分析免费做
我们推文里面提到的各种各样的数据分析环节都是我非常有经验的，比如我在lncRNA的一些基础知识 ，和lncRNA芯片的一般分析流程 介绍过的那些图表，以及下面的目录的分析内容 对我来说是举手之劳，希望可以帮助到你！
转录组数据分析的4个维度认识（数据分析继续免费哦） RNA-seq数据的2个分组差异分析，热图，PCA图，火山图等等
根据感兴趣基因看肝癌免疫微环境的T细胞亚群差异 条形图或者箱线图
查看感兴趣基因的甲基化水平和RNA表达水平(数据分析免费做)相关性 散点图或者箱线图
log与否会改变rpkm形式表达矩阵top的mad基因列表 WGCNA分析免费做
同样的，本次活动我可以帮你免费做一次批量生存分析，但是呢，我也没办法保证结果咋样，有时候数据集就是这样。而且，你需要挑选一下你的阈值哦！
还是老规矩，发送数据分析要求，以及简短的项目描述到我的邮箱 jmzeng1314@163.com
邮件正文最好是加上你是啥时候认识生信技能树的哦，或者其它一些寒暄的话，自我介绍也行。主要是考虑到可能想免费分析数据的朋友很多，所以会根据你的来信，我主观判定一个优先级哦。目前我有20多个愿意长期在我的指导下进行数据探索的学徒，等我的团队扩大到200人，我们应该是可以做到数据分析全部免费，敬请期待哈！