R绘图笔记 | 生存曲线的绘制

104 阅读 0 评论 69 点赞

我是靠谱客的博主哭泣钥匙，这篇文章主要介绍R绘图笔记 | 生存曲线的绘制，现在分享给大家，希望可以做个参考。

参考前文：R绘图笔记 | R语言绘图系统与常见绘图函数及参数

关于绘图，前面介绍了一些：

R绘图笔记 | 一般的散点图绘制

R绘图笔记 | 柱状图绘制

R绘图笔记 | 直方图和核密度估计图的绘制

R绘图笔记 | 二维散点图与统计直方图组合

R绘图笔记 | 散点分布图与柱形分布图

R绘图笔记 | 箱形图的绘制

R绘图笔记 | 小提琴图与漂亮的云雨图绘制

R绘图笔记 | 多数据系列的箱型图与带抖动散点的多数据系列箱型图

R绘图笔记 | 热图绘制

这里介绍生存曲线的绘制，对于生物医学领域来说，文章中很常见的图。

一.数据处理

如果你想获取该数据用于自己练习，下面是获取数据的地址：

https://docs.qq.com/sheet/DV1JxUnpRR3FqbWN3

数据格式是这样的。

数据第A列是病人ID，B到I列临床信息，其他列是病人的标准化后的基因表达数据。

我们只需要左侧红框的临床数据和一列基因表达数据。

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(dplyr)
library(survival)
library(survminer)
data <- read.csv("BioInfoNotesData3.csv",row.names = 1,stringsAsFactors = F)
data <- data[!is.na(data$vital_status),]
data$Time <- ifelse(data$vital_status=="Dead",data$days_to_death,data$days_to_last_followup)
data <- select(data, Time, everything())
surdata <- data[,c(1,4,13)]#选择其中一个基因
colnames(surdata) <- c("Time","Status","Gene")
surdata$Status <- ifelse(surdata$Status=="Dead",0,1)
surdata$Group <- ifelse(surdata$Gene>median(surdata$Gene),"High","Low")
surdata$Time <- surdata$Time/365#最后单位是年
surdata <- surdata[,c(1,2,4)]
head(surdata)

处理后的数据如下：

一列时间，一列生存状态，如果样本的状态是“Dead”，生存分析时候用死亡时间；同样的，如果样本状态是"Alive"，那就用最后随访的时间。Alive为0，Dead为1。一列根据基因表达值高低分为2组。

二.绘图

绘图我们用survival包中的survfit函数，该包是用于生存分析的包。这里简单提一下生存分析里面的参数，参数formula写成Surv(Time,Status) ~ Group 的形式，这里是按照前面的数据写的，你名称不一样可以换，但顺序不能变。当然，该包还可以用于其他模型的分析。

复制代码

1
2
3
4
survfit(formula, data, weights, subset, na.action,  
        stype=1, ctype=1, id, cluster, istate, timefix=TRUE,
        etype, error,  ...)

我们用前面的数据拟合曲线

复制代码

1
2
fit <- survfit(Surv(Time, Status) ~ Group, data = surdata)

拟合的曲线是可以直接用plot来绘图的。

复制代码

1
2
plot(fit, lty = 2:3)

但用plot来绘图需要自己写很多参数，既然有人造好了轮子，我们就直接用就行啦。

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
ggsurvplot(fit,palette = c("#E7B800", "#2E9FDF"),# 自定义调色板
           #conf.int =TRUE,# 增加置信区间
           xlab ="Time in Year", #x轴标题
           urv.median.line = "hv", #增加中位生存时间
           risk.table = TRUE,# 添加风险表
           pval = TRUE, # 添加P值
           legend.title = "",# 设置图例标题，这里设置不显示标题，用空格替代
           legend = c(0.8,0.75), # 指定图例位置
           legend.labs = c("High", "Low")# 指定图例分组标签
           #ggtheme =theme_light() 
)

每一个参数在上面都已经解释了，你可以改改，看看结果图，就知道各参数含义。

另外，还有一个ggsurvplot函数也用来绘制生存曲线。

复制代码

1
2
3
4
5
6
7
8
9
10
ggsurvplot(fit, data = NULL, fun = NULL, color = NULL,
           palette = NULL, linetype = 1, conf.int = FALSE, 
           pval = FALSE, pval.method = FALSE, 
           test.for.trend = FALSE, surv.median.line = "none", 
           risk.table = FALSE, cumevents = FALSE,
           cumcensor = FALSE, tables.height = 0.25, 
           group.by = NULL, facet.by = NULL, add.all = FALSE, 
           combine = FALSE, ggtheme = theme_survminer(), 
           tables.theme = ggtheme, ...

参数解释：

fit # 拟合的生存曲线对象，就是前面拟合的数据
data # 用来拟合生存曲线的数据集，如果没有提供，那么数据将从'fit'对象中提取。
fun # 常用三个字符参数；"event"绘制累计风险概率图 (f(y)=1-y)；"cumhaz"绘制累计风险(HR)，"pct"绘制生存概率(用百分比表示)
color # 设置生存曲线的颜色。
cumevents # 逻辑词，是否添加累积事件表
cumcensor # 逻辑词，是否添加累积删失表
tables.height = 0.25 # 生存曲线图下所有生存表的高度，数值0-1之间
group.by # 包含分组变量名称的字符向量，向量长度≤2
facet.by # 字符向量，指定绘制分面生存曲线的分组变量(应≤2)的名称
ggtheme=theme_survminer() # 设置ggplot2主题，如ggtheme =theme_light()；theme_survminer、theme_cleantable()，也可以是其他包的主题，比如ggpugr包的主题。
tables.theme # 作用于生存表的ggplot2主题名称。
tables.height # 数值(在[0 - 1]中)，指定主要生存图下所有表的高度。
add.all = FALSE # 逻辑词；是否添加总患者生存曲线到主生存图中
palette # 调色板，默认"hue"。可选调色板有"grey", "npg", "aaas", "lancet", "jco", "ucscgb", "uchicago", "simpsons" 和 "rickandmorty"，除此以外你自己可以指定颜色码，比如前面我指定的"#E7B800", "#2E9FDF"。
linetype = 1 # 设置曲线线型。可以按"strata"设置线型；或按数字向量c(1, 2)或按字符向量c("solid", "dashed")设置
conf.int # 逻辑词；默认FASLE；为TRUE则绘制曲线置信区间
pval = FALSE # 逻辑词；为TRUE则将统计检验计算的p值添加到图上；为数字，则直接指定P值大小，如pval = 0.03；为字符串，则添加字符串到图上，如pval = "p-value: 0.031"
pval.method # 逻辑词，是否添加计算p值的统计方法的文本；只有当 pval = TRUE时, 才会在图上添加检验方法文本
test.for.trend # 逻辑词，默认为FALSE；为TRUE则返回趋势p值的检验，趋势检验旨在检验生存曲线的有序差异
surv.median.line # 在中位生存时间点处绘制水平或垂直线的字符向量；可用值有"none"、"hv"、"h"、"v"；其中v绘制垂直线，h绘制水平线。
risk.table = FALSE # 逻辑词，图上是否添加风险表；"absolute" 显示处于风险中的绝对数量；"percentage" 显示处于风险中的百分比数量；"abs_pct" 显示处于风险中的绝对数量和百分比

其实，参数和上面的ggsurvplot函数是一样的。

复制代码

1
2
3
4
5
6
7
8
9
10
ggsurvplot(fit,palette = c("#FF3030", "#0000FF"),
           xlab ="Time in Year", #x轴标题
           urv.median.line = "hv", #增加中位生存时间
           #risk.table = TRUE,# 添加风险表
           pval = TRUE, # 添加P值
           legend.title = "",# 设置图例标题，这里设置不显示标题，用空格替代
           legend = c(0.8,0.75), # 指定图例位置
           legend.labs = c("High", "Low") ,# 指定图例分组标签
           ggtheme =theme_survminer())

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
ggsurvplot(fit,palette = c("#FF3030", "#0000FF"),
           linetype =1,
           xlab ="Time in Year", #x轴标题
           urv.median.line = "hv", #增加中位生存时间
           #risk.table = TRUE,# 添加风险表
           pval = TRUE, # 添加P值
           pval.coord=c(15,0.75),
           legend.title = "",# 设置图例标题，这里设置不显示标题，用空格替代
           legend = c(0.8,0.75), # 指定图例位置
           legend.labs = c("High", "Low") ,# 指定图例分组标签
           ggtheme =theme_test())