r 保留之前曲线_生存曲线居然能够批量绘制了

生信分析第三步：生存曲线批量绘制

各位解螺旋的小伙伴大家好，我是先锋宇，欢迎大家来到每周日的先锋宇专栏，经过前两期推文的学习，很多小伙伴都私信我说从先锋宇助教的专栏很接地气，自己能够开始慢慢处理数据，并且希望先锋宇助教能够继续把这条线走通。听到解螺旋小伙伴积极正向的反馈，小编心理也是非常开心，那么今天咱们继续往下走，我们在前两期推文中完成数据的下载以及差异分析和单因素COX回归，那么今天小编就带大家进行生存曲线的绘制，先锋宇助教还是本着科研效率第一的原则，我们当然不是一次就绘制一条生存曲线，这样与我的风格不符，这次我就来教大家直接批量出图，收获满满的批量生存曲线图。

写在前面

相信很多小伙伴在看文献的时候总是能够看到作者拿一张figure放置多个生存曲线图，不知道大家想过没有，如果作者一张小图一张小图的画，那可能图还没有画完就直接开始拍桌子了。肯定为了提升科研的效率，我们还是希望把这些重复的工作都交给计算机不厌其烦地去做，然后留更多的时间给我们自己享受生活。好啦，言归正传，我们开始今天批量生存曲线绘制的讲解。

代码演示

经过前面两期专栏的处理(没有跟上的同学赶快去前两期专栏看看，打牢基础才能走得更远)，我们现在已经得到了单因素COX回归的结果。

接下来我们筛选p值小于0.05的基因进行保留，这里我们使用dplyr包中的filter函数进行过滤：

uniTab % dplyr::filter(pvalue < 0.05)

接着我们把单因素COX回归有意义的基因再提取出来，因为刚刚得到的是数据库，我们把第一列取出来即可：

unicox_gene 

得到了单因素COX回归有显著统计学意义的基因之后，我们就要开始进行生存曲线绘制即K-M分析，关于K-M分析和COX回归到底有啥区别，大家可以参考风师兄在生信下篇段位3有详细的讲解。如果用我自己的实用的理解那就是，COX分析筛选变量太多的话我们就再加上K-M分析再筛选一次，相当于双重过滤标准；但是如果你COX回归筛选之后就只有几个基因了，那就没有必要再用K-M分析去筛选了，因为筛选了完了你可能就没有基因了。

接下来我们从单因素COX回归的数据框中把pvalue小于0.05的基因提取出来，这里我们使用dplyr包中的select函数，注意这里要记得使用all_off函数把向量放在函数里面，这样才能提取对应的列：

unicoxSig % dplyr::select(1,2,3,all_of(unicox_gene))unicoxSig$futime 

接下来我们进行生存曲线的绘制，首先绘制生存曲线，我们首先需要解压两个强大的包，survival包和survminer包。

library(survival)library(survminer)

首先为了降低难度，我们先来进行一条生存生存曲线的绘制，我们先提取一个基因的表达量：

single_gene % dplyr::select(1,2,3,4)

然后我们构建一个分组文件，根据基因表达量的中位值进行高低两组的划分：

group  median(single_gene[[4]]),               "high", "low")

然后我们计算高低表达两组之间的生存的p值大小：

diff=survdiff(Surv(futime, fustat) ~ group2, data = gene_surv)pValue=1-pchisq(diff$chisq,df=1)if(pValue<0.001){ pValue="p<0.001"}else{ pValue=paste0("p=",sprintf("%.03f",pValue))}

接下来拟合一个生存函数，这里我们使用survfit函数进行拟合：

fit 

最后我们使用ggsurvplot函数来绘制生存曲线，代码参考来自于生信体系课下篇，需要进一步学习的同学可以参看我们的生信体系课，里面有更多丰富的知识等待大家。

ggsurvplot(fit,           data=single_gene,          conf.int=TRUE,          pval=pValue,          pval.size=5,          legend.labs=c("High", "Low"),          xlab="Time (years)",          ylab="Overall survival",          break.time.by = 1,          risk.table.,          palette=c("#d7191c", "#2b83ba"),          risk.table=T,          risk.table.height=.25)

一张可用于文章发表的生存曲线图就绘制好了：

虽然一张绘制好了，但是我们本期的问题还没有解决，我们不仅要一张，我们要很多张。套用一句经典的话就是，只要小孩子才做选择，成年人是全部都要，而且越多越好，图多了我们才有选择的余地。

接下来我们进行批量绘制，批量绘制的原理无非就是循环，而循环就是一列一列循环，然后每一列绘制一个生存曲线。

for(gene in colnames(unicoxSig)[4:ncol(unicoxSig)]){group  median(unicoxSig[[gene]]),               "high", "low")diff=survdiff(Surv(futime, fustat) ~ group, data = unicoxSig)pValue=1-pchisq(diff$chisq,df=1)if(pValue<0.001){ pValue="p<0.001"}else{ pValue=paste0("p=",sprintf("%.03f",pValue))}fit surPlot=ggsurvplot(fit,           data=unicoxSig,          conf.int=TRUE,          pval=pValue,          pval.size=5,          legend.labs=c("High", "Low"),          legend.title=gene,          xlab="Time (years)",          ylab="Overall survival",          break.time.by = 1,          risk.table.,          palette=c("#d7191c", "#2b83ba"),          risk.table=T,          risk.table.height=.25)pdf(file=paste0("surv/",gene,".pdf"), onefile=FALSE, width=6.5, height=5.5)print(surPlot)dev.off()}

写在最后

先锋宇助教每次在跑循环的时候总感觉就是在收获财富，因为每跑一张图，就有可能放到论文里面，然后构成一个完成的figure，希望大家也能和我有同样的感受！希望继续继续关注挑圈联靠公众号，继续关注我的专栏，希望大家都能在这里学有所获，收到满满的干货，好啦，这期的内容就到这里啦，我们下周日再见~

往期传送门：

让生信工作者失业的神器——DrBioRight，真舍不得告诉你！

高效数据清洗！这个R包太强大了！你一定要试试！(文末附赠小彩蛋)

一站式分析R包来了，承包了生信各种分析！太全能了！

搞定这一步，说明你学R有天赋！TCGA数据从下载到差异分析(附代码)

别走，baby，COX森林需要你

欢迎大家关注解螺旋生信频道-挑圈联靠公号~

— END— 撰文丨先锋宇 排版丨四金兄 值班 | 弘  毅 主编丨小雪球