查看原文
其他

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

生信技能树 生信技能树 2022-06-06
前言年前我提出了一个问题:为什么不用TCGA数据库来看感兴趣基因的生存情况就是一篇文章并没有使用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叔致敬)

01获取列表基因的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
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')




03KEGG分析


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芯片的一般分析流程 介绍过的那些图表,以及下面的目录的分析内容 对我来说是举手之劳,希望可以帮助到你!

同样的,本次活动我可以帮你免费做一次批量生存分析,但是呢,我也没办法保证结果咋样,有时候数据集就是这样。而且,你需要挑选一下你的阈值哦!

还是老规矩,发送数据分析要求,以及简短的项目描述到我的邮箱 jmzeng1314@163.com

邮件正文最好是加上你是啥时候认识生信技能树的哦,或者其它一些寒暄的话,自我介绍也行。主要是考虑到可能想免费分析数据的朋友很多,所以会根据你的来信,我主观判定一个优先级哦。目前我有20多个愿意长期在我的指导下进行数据探索的学徒,等我的团队扩大到200人,我们应该是可以做到数据分析全部免费,敬请期待哈!

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存