我以前写过一个帖子
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
