使用R语言复现Nature Cancer论文中的生存分析

小明的数据分析笔记本

论文

Concurrent delivery of immune checkpoint blockade modulates T cell dynamics to enhance neoantigen vaccine-generated antitumor immunity

论文中提供了大部分的作图的原始数据，今天的论文利用论文中提供的原始数据复现一下论文的fig6c

生存分析需要用到的数据，某种干预后的时间，经过这个时间后的状态一般是存活或者死亡 存活是代码0，死亡是代码1，通常还有一个分组，可以比较干预措施在两组之间的效果是否有差别

读取数据

library(readxl)
library(tidyverse)

read_excel("2024.data/20240629/fig6c.xlsx") %>% 
  select(Time,status,CAST_score,vital_status) -> fig6c.dat

这个数据集里有很多列，指保留可能用到的列

生存分析R语言里有现成的R包 survival

作图也有专门的R包ggsurvfit，这个R包里提供了做生存分析的函数 survfit2()，还提供了函数 tidy__survfit()可以直接把生存分析的结果转化成数据框用于作图

生存分析并组图的代码

survfit2(Surv(Time, status) ~ `CAST_score`, data = fig6c.dat) %>% 
  tidy_survfit() %>% 
  ggplot(aes(x=time,y=estimate))+
  geom_point(aes(color=strata))+
  geom_line(aes(group=strata,color=strata),
            linewidth=2)+
  theme_bw(base_size = 15)+
  theme(panel.grid = element_blank(),
        legend.position = c(0.8,0.8))+
  scale_color_manual(values = c("#56b4e8","#d55e00"),
                     name=NULL)+
  annotate(geom = "text",x=1000,y=0.2,
           label="italic(P)==0.007~(n==214)",
           parse=TRUE,size=8)+
  labs(x="Time (d)",y="Overall survival",
       title = "High CD8 TIL (SKCM)")+
  scale_y_continuous(limits = c(0,1),
                     breaks = seq(0,1,by=0.2))+
  scale_x_continuous(breaks = seq(0,10000,by=2000))+
  theme(axis.text.x = element_text(angle=60,hjust=1,vjust=1))

计算p值的代码

survdiff(Surv(Time, status) ~ `CAST_score`, data = fig6c.dat)

这里做的是卡方检验，如何得出论文中图上呈现的FDR值暂时想不到办法了。

参考链接

• 1、https://bioconnector.github.io/workshops/r-survival.html
• 2、https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html
• 3、https://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/R/R7_LogisticRegression-Survival/R7_LogisticRegression-Survival4.html
• 4、https://github.com/pharmaverse/ggsurvfit/blob/main/R/ggsurvfit.R

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！

微信扫一扫分享文章