8秒完成2万个基因的生存分析，人人都可以！

我以前写过一个帖子
TCGA真实数据下的批量生存分析

心里很高兴，因为这是我学习生信后做的第一件像样的事，解决了我心里多年的麻烦。当时，正常运行20000个基因要花费50分钟。但是，今天，我10s钟就实现了。

事情的经过是这样的。
首先我们加载生存数据,也可以通过上次那个帖子来准备

rm(list = ls())
library(survival)
load("survival_data.Rdata")
rt <- survival_data
test <- rt[1:10,1:10]

这个数据是个数据框，有209行，有19452列，截取一部分来看看

前面两列分别是生存状态，生存时间，后面第3列开始全是基因的表达量。
这个格式很重要，无论是生存分析，还是模型构建，都是这个格式。
这个数据框的名字现在叫做rt(单存为了简单而简单)

for循环批量生存分析

一到批量，我脑子里面跳出的是for循环，因为，我觉得他无论如何都会完成我的任务。代码是这个样子的。

## 创建空的数据框
res2 <- data.frame()
## 获取基因列表
genes <- colnames(rt)[-c(1:2)]
## for循环开始
for (i in 1:length(genes)) {
  print(i)
  # 中位数分组
  group = ifelse(rt[,genes[i]]>median(rt[,genes[i]]),"high","low")
  surv =as.formula(paste('Surv(futime, fustat)~', "group"))
  data = cbind(rt[,1:2],group)
  #生存分析求差异
  x = survdiff(surv, data = data)
  # 获取p值
  pValue=1-pchisq(x$chisq,df=1) 
  #第一列基因名称
  res2[i,1] = genes[i]
  #第二列是p值
  res2[i,2] = pValue
}

运行到1374个的时候报错了，我们说，报错的时候不要怕，要去读。

意思是说，这个基因，无法分为两组，也就是说，这个基因的表达量都一样，这在真实状态下是不可能的，应该是本身表达量很低，标准化后变成一样的了。

所以，我们需要添加语句告诉他，碰到这样的基因，请你跳过去。这个事情是if做的。

if(length(table(group))==1) next

添加if判断语句，再用system.time函数计算速度，正是运行代码如下

res2 <- data.frame()
genes <- colnames(rt)[-c(1:2)]
system.time(for (i in 1:length(genes)) {
  print(i)
  group = ifelse(rt[,genes[i]]>median(rt[,genes[i]]),"high","low")
  if(length(table(group))==1) next
  surv =as.formula(paste('Surv(futime, fustat)~', "group"))
  data = cbind(rt[,1:2],group)
  x = survdiff(surv, data = data)
  pValue=1-pchisq(x$chisq,df=1) 
  res2[i,1] = genes[i]
  res2[i,2] = pValue
})
names(res2) <- c("ID","pValue_log")

最终用时50秒，比以前要快很多，当时是用了50分钟。

lapply批量生存分析

for循环能达到这个速度，那么lapply速度肯定会更快，我们把for循环改成function，配合lapply来试一下。这里面有一个限速点，是，刚才的if判断next失效了，因为现在是函数，所以我们把next改成return(NULL)就可以了。

完整的代码是这个样子的。

genes <- colnames(rt)[-c(1:2)]
system.time(res3 <- lapply(1:length(genes), function(i){
  group = ifelse(rt[,genes[i]]>median(rt[,genes[i]]),"high","low")
  if(length(table(group))==1) return(NULL)
  surv =as.formula(paste('Surv(futime, fustat)~', "group"))
  data = cbind(rt[,1:2],group)
  x = survdiff(surv, data = data)
  pValue=1-pchisq(x$chisq,df=1) 
  return(c(genes[i],pValue))
}))

速度变成了34秒，有进步。

future.apply批量生存分析

既然已经有了函数，那么还可以考虑一下，并行化处理，foreach等包可以实现，不过我觉得太麻烦了，有那个时间，还不如单存使用lapply呢。有没有什么方法，几乎不需要改代码，实现并行化操作呢？

有！future.apply就可以，对于apply家族的成员，apply，lapply，mapply, tapply等，只要把他们写成future_apply, future_lapply，future_mapply, future_tapply就可以了。实在是懒人福音。
使用起来也很方便，首先加载R包

library(future.apply)

其次，告诉R语言我要并行化

plan(multiprocess)

最后，原来的lapply改成future_lapply

运行一下，真的是8s以内！

lapply返回的是list，需要用do.call来转换

res4 <- data.frame(do.call(rbind,res4))
names(res4 ) <- c("ID","pValue_log")

按照p从小到大排序，选取出p值小于0.05的基因,有1453个。

library(dplyr)
res5 <- res4 %>%
        filter(pValue_log < 0.05)%>%
        arrange(pValue_log)

选取第二个作图了解一下：

index <- res5$ID[2]
group = ifelse(rt[,index]>median(rt[,index]),"high","low")
surv =as.formula(paste('Surv(futime, fustat)~', "group"))
data = cbind(rt[,1:2],group)
my.surv <- Surv(rt$futime, rt$fustat)
fit <- survfit(my.surv ~ group)
library(survminer)
ggsurvplot(fit, data = data)

这只是其中一个。。

总结一下：

• 批量的事情要想到写for循环
• 如果特殊情况要处理，要用if语句
• 偷懒并行化请用future.apply

更多帖子，关注“果子学生信”微信公众号。