使用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))

image.png
计算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
欢迎大家关注我的公众号
小明的数据分析笔记本