生信分析、meta分析、数据挖掘
TCGA、GEO、SEER、Oncomine

甜过初恋!这次是真的批量做TCGA的生存分析

大家好,今天是12月4号,农历10月17,额,好像今天除了要开组会以外,也不是个什么特别的日子。

在刚刚进入生信领域的时候,我想做的事情就是三个, 

第一知道任何我想研究的基因在组织中的表达情况,

第二,我选的基因对肿瘤的生存有无影响, 

第三这个基因可能的作用是什么?

这是来自临床医生的视角,研究疾病,最终希望能够服务临床,临床离不开诊断和治疗,假设一个基因的表达对肿瘤的预后有影响,他很可能就是我的盘中餐。

有一大堆网页工具可以实现生存分析,但是你看看jimmy已经写的帖子

都可以批量做生存分析了,还要网页工具干嘛?

一拳把人打翻在地,扶都扶不起来。 但是没有办法他是群主。即使会随时被朝阳群众扭送到派出所,他依然是群主,他的排版有问题,但是架不住他的内容有深度,群主喜闻乐见。

那么问题来了,上面的问题也可以反过来问,既然别人已经准备好了网页工具给你用,你还写代码做什么?!

四个字,批量,自由

大气一点就是高端定制。

就这么简单。 那开始我们的表演:

今天我们要对TCGA里面的任意基因做生存分析,最关键的我们要批量做生存分析,然后选取生存差异最显著的基因。

首先安装需要的包:

  1. # Load the bioconductor installer.

  2. source("https://bioconductor.org/biocLite.R")

  3. options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

  4. # Install the main RTCGA package

  5. biocLite("RTCGA")

  6. # Install the clinical and mRNA gene expression data packages

  7. biocLite("RTCGA.clinical")

  8. biocLite("RTCGA.mRNA")

然后是加载包

  1. library(RTCGA)

  2. #了解数据

  3. infoTCGA <- infoTCGA() #这个命令会返回一个数据框,可以知道有哪些数据可被下载

  4. #获得临床数据:

  5. # Create the clinical data

  6. library(RTCGA.clinical)

  7. clin <- survivalTCGA(BRCA.clinical) #到这里临床部分的信息已经获得啦

得到数据后我们先看一下他是什么结构

  1. class()

[1] “data.frame”

再看一下前面几行数据

  1. head(clin)

times bcrpatientbarcode patient.vital_status 1 3767 TCGA-3C-AAAU 0 2 3801 TCGA-3C-AALI 0 3 1228 TCGA-3C-AALJ 0 4 1217 TCGA-3C-AALK 0 5 158 TCGA-4H-AAAK 0 6 1477 TCGA-5L-AAT0 0

简单说就是三列,TCAGid,生存时间,和发生的事件

获得gene表达数据:

  1. library(RTCGA.mRNA) #加载数据包

  2. class(BRCA.mRNA)  #查看数据类型发现是个数据框

  3. dim(BRCA.mRNA)  #看一下数据维度发现有590个样本,17815个基因

  4. BRCA.mRNA[1:5, 1:5] #看一下部分数据样子

  1.      bcr_patient_barcode    ELMO2  CREB3L1    RPS11  PNMA1

1 TCGA-A1-A0SD-01A-11R-A115-07 0.5070833 1.43450 0.765000 0.52600 2 TCGA-A1-A0SE-01A-11R-A084-07 0.1814167 0.89075 0.716000 0.13175 3 TCGA-A1-A0SH-01A-11R-A084-07 0.4615000 2.25925 0.417125 0.32500 4 TCGA-A1-A0SJ-01A-11R-A084-07 0.8770000 0.43775 0.115000 0.75775 5 TCGA-A1-A0SK-01A-12R-A084-07 1.4123333 -0.63725 0.492875 0.94325

好了,行是样本,列为gene 下面我们挑选几个基因的表达数据出来,融合到生存的数据上去,因为表达量数据中TCGA的id号要长一点 所以在融合前,需要先裁剪一下,为了避免产生过多的中间变量,我们使用管道符号%>%,他的作用是 把前一个计算得到结果,作为第二个函数的参数,示例如下:

  1. library(dplyr)

  2. exprSet <- BRCA.mRNA %>%

  3.  # then make it a tibble (nice printing while debugging)

  4.  as_tibble() %>%

  5.  # then get just a few genes,这里是测试用

  6.  select(bcr_patient_barcode, PAX8, GATA3, ESR1) %>%

  7.  # then trim the barcode (see head(clin), and substr)

  8.  mutate(bcr_patient_barcode = substr(bcr_patient_barcode, 1, 12)) %>%

  9.  # then join back to clinical data

  10.  inner_join(clin, by="bcr_patient_barcode")

看一下数据,发现就是表达量加上time,event两列,所以记住这规律,最终批量的时候需要减掉这两列


mark


开始做生存分析

  1. library(survival)

  2. library(survminer)

对需要做生存分析的样本分组,把连续变量变成分类变量,这里选择测试的基因是GATA3

  1. group <- ifelse(esprSet$GATA3>median(esprSet$GATA3),'high','low')

构建生存对象,并且进行数据处理,这里是两步,我只是融合在了一起

  1. sfit <- survfit(Surv(times, patient.vital_status)~group, data=esprSet)

  2. sfit

  3. summary(sfit)

直接绘图

  1. ggsurvplot(sfit, conf.int=F, pval=TRUE)


mark



单个基因绘制成功,我看一下能不能小规模循环操作 首先得到得到生存对象my.surv

  1. my.surv <- Surv(esprSet$times, esprSet$patient.vital_status)

使用apply循环,对数据的2到4列进行操作,实际上就是PAX8, GATA3, ESR1这三个基因, 这里面使用了survdiff,用来比较差异大小,获得p值

  1. log_rank_p <- apply(esprSet[,2:4], 2, function(values1){

  2.  group=ifelse(values1>median(values1),'high','low')

  3.  kmfit2 <- survfit(my.surv~group,data=esprSet)

  4.  #plot(kmfit2)

  5.  data.survdiff=survdiff(my.surv~group)

  6.  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)

  7. })

运行完了之后返回的找出小于0.05的P.val,在这个例子里面因为没有基因的p.val小于0.05,所以筛选不出来

  1. log_rank_p <- log_rank_p[log_rank_p<0.05]

筛选后排序,并获得基因名 genediff <- as.data.frame(sort(logrank_p))


好了生存数据已经获取,

表达数据小规模试验可行,

批量操作也能实现,

下面就进入实战环节,前戏实在太长,但是你要相信你只要不睡着,这些都是值得的

获得完整的表达量数据

  1. library(dplyr)

  2. esprSet <- BRCA.mRNA %>%

  3.  # then make it a tibble (nice printing while debugging)

  4.  as_tibble() %>%

  5.  # then trim the barcode (see head(clin), and ?substr)

  6.  mutate(bcr_patient_barcode = substr(bcr_patient_barcode, 1, 12)) %>%

  7.  # then join back to clinical data

  8.  inner_join(clin, by="bcr_patient_barcode")

构建生存对象my.surv

  1. library(survival)

  2. my.surv <- Surv(esprSet$times, esprSet$patient.vital_status)

在进行下一步之前,我居然突发奇想,我想看一看,这个表达数据里面哪些基因的NA值最多,有没有NA值多过样本数量一般的基因呢? 先构建了一个函数,他对数据的列起作用,统计NA值的个数,最终返回成一个数据框

  1. rem <- function(x){

  2.  r <-as.numeric(apply(x,2,function(i) sum(is.na(i))))

  3.  return(data.frame(geneName=names(x)[which(r > 0)],na_num=r[which(r > 0)]))

  4. }

然后对表达量数据进行统计

  1. na_count <- rem(esprSet)

最终发现NA最多的基因是LCE1B,有17个,所以数据不需要特殊处理啦

  1. na_count <- dplyr::arrange(na_count,desc(na_num))

那就开始批量运算了,一开始就用apply,发现大概需要运行50分钟以上,所以尝试使用并行化处理 R语言里面的并行化有个专门的项目就是给apply的,使用起来也是很方便

  1. #尝试使用并行运算

  2. library(parallel)

  3. #detectCores()检查当前电脑可用核数

  4. cl.cores <- detectCores()

  5. #makeCluster(cl.cores)使用刚才检测的核并行运算,我的服务器是28核56线程,我就用50吧

  6. cl <- makeCluster(50)

  7. #这是坑,parApply里面用到的函数以及变量都需要申明,不声明就必须用模块

  8. clusterExport(cl,c("esprSet","my.surv"))

  9. #length(names(esprSet))-2,为什么减去2,因为之前小规模测试时,我们知道最后两个是time和event,不是表达量

  10. #数据从25开始,原因是从2开始会报错,暂时无法解决,还有要注意是parApply,A要大写的

  11. log_rank_p <- parApply(cl,esprSet[,25:length(names(esprSet))-2],2,function(values){

  12.  group=ifelse(values>median(na.omit(values)),'high','low')

  13.  kmfit2 <- survival::survfit(my.surv~group,data=esprSet)

  14.  #plot(kmfit2)

  15.  data.survdiff=survival::survdiff(my.surv~group)

  16.  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)

  17. })

这个运行的时间可能就是2分钟 终止并行化

  1. stopCluster(cl)

找出小于0.05的P.val

  1. log_rank_p <- log_rank_p[log_rank_p<0.05]

筛选后排序,并获得基因名

  1. gene_diff <- as.data.frame(sort(log_rank_p))

最终得到的gene是2153个,这个数据是我之前留下的,本次写贴时要带孩子,在家没法运行运算。 数据保存可以这样:

  1. save(gene_diff,file = "gene_df.Rda")

如果想用的时候就这样:

  1. load(file = "gene_df.Rda")

没有必要先写成txt格式,然后要用的时候再读取进来,直接保存成R语言的格式即可,

你能想象每次用完word保存成图片,下次再使用时用OCR识别图片变成文字再编辑的状态么?

既然到了这一步,我们随便选取一个基因来作图,闭着眼睛都知道,都是有差异的

  1. library(survminer)

  2. group <- ifelse(esprSet$LRRC8D>median(esprSet$LRRC8D),'high','low')

  3. sfit <- survfit(Surv(times, patient.vital_status)~group, data=esprSet)

  4. ggsurvplot(sfit, conf.int=FALSE, pval=TRUE)


mark

好吧效果很不错嘛

这时候把癌和癌旁的数据作差异分析,得到的基因与今天获得的基因取交集,就可以获得又差异表达,又对生存有影响的基因了。

今天我们是用的别人已经下载好的数据,明天我们来尝试自己下载并且清理数据,而且只用一种语言,就是R语言

要知道,今天往后的ceRNA网络构建,单基因GSEA都是基于这些表达量数据的。

今天很美好,明天更有用,但是大部分人都。。。。马云说的。


赞(0) 打赏
未经允许不得转载:医学SCI科研之家 » 甜过初恋!这次是真的批量做TCGA的生存分析
分享到: 更多 (0)

评论 抢沙发

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址

meta分析、生信分析

meta、生信交流群综合科研交流群