R语言——绘制生命曲线图（细胞因子IL5）

1. R语言的官方环境
   R Base
   官网下载：
   https://www.r-project.org/
   包含R的核心解释器和基础功能，适合纯命令行操作。
   RStudio
   官网下载：
   https://www.rstudio.com/products/rstudio/download/

说明

初学者注意
（1）代码中所有符号不能是中文
（2）#表示后面的内容是注释，不生效
（3）括号是成对出现的

关于NA值

在 R 语言中，NA（Not Available）是用于表示缺失值（Missing Value）的特殊值。它表示某个数据点不存在、不可用或未被记录。

1. NA 在计算中的行为
   大部分数学函数（如 sum(), mean()）默认会因 NA 返回 NA
   需使用 na.rm = TRUE 忽略 NA 计算

2. 如何检查NA值：
   用 is.na()

3. 如何处理NA值：
   （1）删除 NA
   na.omit()：删除所有含 NA 的行
   complete.cases() + 索引：筛选完整数据
   （2）替换 NA
   均值/中位数填充（适用于数值型数据）
   众数填充（适用于分类数据）
   插值法（时间序列数据）
   （3）计算时忽略 NA
   使用 na.rm = TRUE 忽略 NA 计算

流程

（1）加载需要的依赖
（2）读取文件
（3）预处理文件：保障所需数据的完整性、格式规范；
（4）创建对象，拟合生存曲线的数据
（5）绘制曲线

代码

加载包

library(survival)
library(survminer)
library(dplyr)
library(readxl)

如果提示没有包，就需要下载，例如：不存在叫‘survival’这个名称的程序包

install.packages("survival")

读取Excel文件

data <- read_excel("G:/术前.xlsx", sheet = "术前")

"G:/术前.xlsx"是文件名，"术前"是工作区名称

清理数据

重命名列名

# 重命名列名
colnames(data) <- c("IL5", "Status", "OS_months")

c表示列名,这里没有写列序号，是按顺序修改前三列，如果要指定修改某列列名，示例如下：

# 方法1：直接通过列索引修改特定列名
colnames(df)[3] <- "name"
colnames(df)[7] <- "data"
colnames(df)[11] <- "time"
# 方法2：使用向量一次性修改多个列名（更简洁）
colnames(df)[c(3,7,11)] <- c("name", "data", "time")

处理IL-5中的"<"符号 - 替换为检测下限的一半

首先找出所有包含"<"的值

data$IL5_clean <- ifelse(grepl("<", data$IL5), as.numeric(gsub("[^0-9.]", "", data$IL5))/2,as.numeric(data$IL5))

函数	说明	特点
ifelse(test, yes, no)	test：逻辑测试条件，yes：当test为TRUE时返回的值，no：当test为FALSE时返回的值	可以同时处理整个向量/列，会保持输入数据的结构和属性
grepl(“<”, data$IL5	在IL5中搜索"<"符号	默认区分大小写，可设置ignore.case=TRUE；多列中查找c(“12”, “<5”, “3<2”, “NA”)
as.numeric(data$IL5)	将IL5整列的数据强制转换为数值型	无法转换时会生成NA并给出警告
gsub(“[^0-9.]”, “”, data$IL5)	把IL5中所有不是数字和小数点的字符串替换为空	应用场景：非常适合清洗实验数据中的检测限值（如"<0.05"），可用于提取混杂在文本中的数值（如"23.5mg/L"），常用于处理实验室数据或医学检测报告中的数值
data$IL5_clean	如果IL5_clean列不存在，会创建这个新列	详细说明请看补充部分：data$new_column 创建新列/修改现有列

注意：
这个清理过程，可能会有 某些值无法被正确转换为数值的情况：比如数据中存在空字符串、非数字字符或其他格式问题
解决方法请看补充部分：更健壮的数据清理方法

检查缺失

有些行的数据不合规，通过下述代码可以进行统计不合规的数据有多少

sum(is.na(data$IL5_clean))
sum(is.na(data$Status))
sum(is.na(data$OS_months))

sum(is.na(data$ IL5_clean))：统计IL5_clean列的缺失值数量
sum(is.na(data$ Status))：统计Status列的缺失值数量
sum(is.na(data$ OS_months))：统计OS_months列的缺失值数量
输出示例：
[1] 5 # IL5_clean有5个缺失值
[1] 2 # Status有2个缺失值
[1] 3 # OS_months有3个缺失值

移除缺失值

简单粗暴的方式，把不合规的数据移除（移除整行）

data_clean <- data %>% filter(!is.na(IL5_clean) & !is.na(Status) & !is.na(OS_months))

创建一个新的数据框data_clean，只保留三个关键列都非缺失值的完整记录。

data_clean <- data：将数据框data复制一份到data_clean
管道操作符%>% 作用：将左侧对象传递给右侧函数的第一个参数
filter(.data, …, .preserve = FALSE)

.data：要筛选的数据框
…：逻辑条件表达式
.preserve：是否保留分组结构（高级用法）

意义

1. 数据质量评估：
   先统计缺失值了解数据完整度
   判断是否需要插补或可以直接删除

2. 分析准备：
   许多统计方法要求完整数据（如生存分析）
   确保后续分析基于相同的一组观测值

IL5_clean：生物标志物测量值
Status：患者状态（如生死）
OS_months：总生存期
三者都是关键分析变量，必须完整

根据IL-5中位数将患者分为高低两组

# 根据IL-5中位数将患者分为高低两组
median_il5 <- median(data_clean$IL5_clean)
data_clean$IL5_group <- ifelse(data_clean$IL5_clean > median_il5, "High", "Low")

median() 计算中位数，默认情况下，如果向量包含NA值，函数会返回NA，除非设置na.rm = TRUE
median_il5 <- 将计算得到的中位数赋值给新变量median_il5，供后续使用。
data_clean$IL5_group <- 将ifelse的结果赋值给数据框中新创建的列IL5_group。

结果呈现：

创建生存对象

surv_obj <- Surv(time = data_clean$OS_months, event = data_clean$Status)

Surv() 创建生存数据对象，这个对象将作为后续生存分析函数的输入。

注意1：
Surv() 函数要求时间变量必须是数值型。如果检查出变量类型不符合，会返回NA值
解决方案请看补充部分：类型转换（转换为数值型）
注意2：
event 参数必须是一个逻辑值（TRUE/FALSE）或数值型（0/1），其中：
0 或 FALSE 表示 删失（censored）（患者存活或失访）
1 或 TRUE 表示 事件发生（死亡）
解决方案请看补充部分：类型转换（转换为数值型0/1）

拟合生存曲线

fit <- survfit(surv_obj ~ IL5_group, data = data_clean)

survfit() 拟合生存曲线，默认使用Kaplan-Meier估计法。

参数解释
surv_obj ~ IL5_group:
R公式(formula)，指定了生存分析的结构
左边surv_obj是我们之前用Surv()创建的生存对象
右边IL5_group是分组变量(基于之前IL5中位数分组的High/Low组)
data = data_clean:
指定了变量所在的数据框
这样可以直接使用列名(IL5_group)而不需要用data_clean$IL5_group

注意：关于如何保存拟合后的生存曲线对象，请看补充部分：保存 / 加载 拟合后的生存曲线对象

绘制生存曲线

# 使用ggsurvplot函数绘制Kaplan-Meier生存曲线图
surv_plot<-ggsurvplot(# 必需参数：生存分析结果对象，通常由survfit()函数生成fit, # 指定用于生存分析的数据集data = data_clean,# 设置y轴显示范围，为p值标签留出空间#ylim = c(0.9, 1.0),  # 更合理的范围（留出标签空间）# 逻辑值，是否在图上显示log-rank检验的p值(默认显示在右上角)pval = TRUE, # 逻辑值，是否显示p值的方法标签(如"Log-rank test")pval.method = TRUE,# 设置p值标签的位置
pval.method.coord = c(1, 1.0), # 方法标签放在 (x=1, y=1.0)pval.coord = c(2, 1.0),       # 将 P 值放在 (x=2, y=1.0)# 是否显示生存曲线的置信区间带conf.int = FALSE,# 设置置信区间带的透明度（0-1，越小越透明）# conf.int.alpha = 0.2,          # 20%透明度# 逻辑值，是否在图形下方添加风险表(显示各时间点的风险人数)risk.table = TRUE,# 设置风险表高度（占整个图形高度的比例）risk.table.height = 0.25,      # 占25%的高度# 自定义图例标签：# 使用paste()动态创建分组标签，基于IL-5的中位数分组# 第一组标签："Low IL-5 (≤中位数值)"# 第二组标签："High IL-5 (>中位数值)"# round(median_il5, 2)将中位数保留两位小数legend.labs = c(paste("Low IL-5 (≤", round(median_il5, 2)), paste("High IL-5 (>", round(median_il5, 2))),# 设置图例的标题为"IL-5 Group"legend.title = "IL-5 Group",# 设置x轴标签为"Time (Months)"，表示时间以月为单位xlab = "Time (Months)",# 设置y轴标签为"Survival Probability"，表示生存概率ylab = "Survival Probability",# 设置图形的主标题title = "Kaplan-Meier Survival Curve by IL-5 Level",# 设置两条生存曲线的颜色：# 第一组颜色为十六进制码#E7B800(红色)# 第二组颜色为十六进制码#2E9FDF(蓝色)palette = c("#FF6A6A", "#2E9FDF"),# 设置图形主题为ggplot2的简约主题(去除多余背景元素)ggtheme = theme_minimal()
)

RGB颜色对照表

如何保存绘制出来的曲线，请看补充部分：保存绘制的生存曲线（如图片或PDF）

补充

data$new_column 创建新列/修改现有列

一、列操作 (使用$和[ ])

1. 单列操作

# 使用$选择单列（返回向量）
data$column_name# 使用[ ]选择单列（返回数据框）
data["column_name"]
data[, "column_name"]  # 逗号在前表示选择所有行

2. 多列操作

# 选择多列（返回数据框）
data[c("col1", "col2", "col3")]
data[, c("col1", "col2")]# 使用列索引
data[, 1:3]  # 选择前3列

3. 整列操作

# 创建/替换整列
data$new_column <- values_vector  # 长度需匹配行数# 删除列
data$column_name <- NULL
data[, "column_name"] <- NULL
data <- subset(data, select = -column_name)

二、行操作 (使用[ , ])

1. 单行操作

# 选择单行（返回数据框）
data[1, ]  # 第一行
data["row_name", ]  # 如果有行名# 修改单行
data[1, ] <- c(value1, value2, ...)

2. 多行操作

# 选择多行
data[1:5, ]  # 1-5行
data[c(1,3,5), ]  # 第1,3,5行# 条件选择
data[data$age > 30, ]  # age大于30的行
subset(data, age > 30)  # 等效方法

3. 整行操作

# 添加新行
data <- rbind(data, new_row)# 删除行
data <- data[-c(1,3), ]  # 删除第1和3行
data <- data[data$age >= 18, ]  # 删除age<18的行

三、混合操作

1. 选择特定行列

# 选择第1-3行，第2-4列
data[1:3, 2:4]# 条件选择特定列
data[data$age > 30, c("name", "age")]

2. 修改特定区域

# 修改第2-3行的"score"列
data[2:3, "score"] <- c(90, 85)# 使用逻辑条件修改
data[data$gender == "M", "salary"] <- data[data$gender == "M", "salary"] * 1.1

更健壮的数据清理方法

data$IL5_clean <- sapply(data$IL5, function(x) {if (is.na(x) || x == "") {return(NA)  # 如果是NA或空字符串，返回NA} else if (grepl("<", x)) {# 提取数值并除以2num <- as.numeric(gsub("[^0-9.]", "", x))if (is.na(num)) {return(NA)  # 如果提取失败，返回NA} else {return(num / 2)}} else {# 尝试直接转换为数值num <- as.numeric(x)return(ifelse(is.na(num), NA, num))}
})

额外建议
在转换前先检查数据的唯一值：

unique(data$IL5)

（1）如果数据中有很多不规则的值，可以考虑手动处理特殊情况
（2）对于无法自动转换的值，可以单独处理或考虑是否应该排除这些观测
（3）这样修改后，应该能更稳健地处理数据转换问题，并清楚地知道哪些值无法被转换。

处理完后检查结果：

# 检查转换后的结果
sum(is.na(data$IL5_clean))  # 查看有多少NA值
which(is.na(data$IL5_clean))  # 查看哪些行有NA值
data[is.na(data$IL5_clean), ]  # 查看这些行的原始数据

类型转换（转换为数值型）

> surv_obj <- Surv(time = data_clean$OS_months, event = data_clean$Status)
错误于Surv(time = data_clean$OS_months, event = data_clean$Status): Time variable is not numeric

data_clean$OS_months 不是数值型（numeric）变量，而 Surv() 函数要求时间变量必须是数值型。我们需要先检查并转换这个变量的类型。

转换方法：

1. 检查变量类型

# 查看变量类型
str(data_clean$OS_months)
class(data_clean$OS_months)

2. 转换变量类型
   如果 OS_months 是字符型（character）或因子型（factor），需要转换为数值型：

# 转换为数值型
data_clean$OS_months <- as.numeric(as.character(data_clean$OS_months))# 再次检查转换结果
sum(is.na(data_clean$OS_months))  # 查看是否有转换失败的NA值
which(is.na(data_clean$OS_months))  # 查看哪些行转换失败

3. 处理转换失败
   如果转换产生了 NA 值，可能是因为原始数据中有非数字字符（如"／"或其他符号）：

# 查看原始数据中的特殊值
unique(data$OS_months)[!is.na(unique(data$OS_months)) & is.na(as.numeric(as.character(unique(data$OS_months))))]# 清理特殊字符（例如将"／"替换为NA）
data_clean$OS_months <- ifelse(data_clean$OS_months %in% c("／", "/", " "), NA, data_clean$OS_months)
data_clean$OS_months <- as.numeric(as.character(data_clean$OS_months))# 移除NA值
data_clean <- data_clean[!is.na(data_clean$OS_months), ]

然后重新转换

# 确保Status是数值型（0/1）
data_clean$Status <- as.numeric(as.character(data_clean$Status))# 再次检查转换结果
sum(is.na(data_clean$OS_months)) 
which(is.na(data_clean$OS_months))  

完成后，查看的结果如下，可以看到已经是我们需要的整数类型

类型转换（转换为数值型0/1）

> surv_obj <- Surv(time = data_clean$OS_months, event = data_clean$Status)
错误于Surv(time = data_clean$OS_months, event = data_clean$Status): Invalid status value, must be logical or numeric
>

解决方法：

1. 检查 Status 的当前值和类型

# 查看 Status 列的唯一值和类型
unique(data_clean$Status)
class(data_clean$Status)

2. 转换 Status 为合法的 0/1 格式

(1) 如果 Status 是字符型（如 “0”, “1”）

data_clean$Status <- as.numeric(as.character(data_clean$Status))

(2) 如果 Status 是字符串（如 “生存”, “死亡”）

data_clean$Status <- ifelse(data_clean$Status == "死亡", 1, 0)

(3) 如果 Status 包含其他值（如 2, NA）

# 只保留 0 和 1，其余设为 NA 并移除
data_clean$Status <- ifelse(data_clean$Status %in% c(0, 1), data_clean$Status, NA)
data_clean <- data_clean[!is.na(data_clean$Status), ]

3. 重新查看

 # 查看 Status 列的唯一值和类型
unique(data_clean$Status)
class(data_clean$Status)

保存 / 加载 拟合后的生存曲线对象

1. 说明：
   通过survfit()拟合的生存曲线对象，包含以下关键信息：
   生存概率：每个时间点的生存率（surv）。
   时间点：事件发生的具体时间（time）。
   风险集：每个时间点处于风险中的患者数（n.risk）。
   事件数：每个时间点发生事件的数量（n.event）。
   分组信息（如果存在，如本例中的IL5_group）。

在R中，保存拟合的生存曲线对象（如fit）后，可以随时重新加载并使用它进行绘图或分析

1. 保存生存曲线对象
   使用save()或saveRDS()函数将对象保存到本地文件：
   方法1：使用 save()（保存为RData格式）

# 保存单个对象
save(fit, file = "survival_fit.RData")# 保存多个对象（如同时保存surv_obj和fit）
save(surv_obj, fit, file = "survival_objects.RData")

优点：可保存多个对象到一个文件。
文件格式：二进制（.RData或.rda）。

方法2：使用 saveRDS()（保存为单个对象）

saveRDS(fit, file = "survival_fit.rds")

优点：更灵活，加载时需指定新变量名（见下文）。
文件格式：二进制（.rds）

2. 加载保存的对象
   根据保存方式选择对应的加载函数：

加载 save() 保存的文件

load("survival_fit.RData")  # 直接恢复原对象名（如fit）

加载后，对象会以原名（如fit）自动出现在环境中。

加载 saveRDS() 保存的文件

fit_loaded <- readRDS("survival_fit.rds")  # 需赋值给新变量名

可自由命名加载后的对象（如fit_loaded）。

3. 使用加载的对象绘制生存曲线
   加载后的对象和原始对象完全一致，可直接用于绘图：

使用 survminer 包绘图

library(survminer)
ggsurvplot(fit_loaded,           # 加载后的对象data = data_clean,    # 需确保数据存在pval = TRUE, risk.table = TRUE,conf.int = TRUE,palette = "jco")      # 自定义颜色

使用基础R绘图

plot(fit_loaded, col = c("red", "blue"),xlab = "Time (months)", ylab = "Survival Probability")
legend("topright", legend = levels(data_clean$IL5_group), col = c("red", "blue"), lty = 1)

4. 注意事项
   数据一致性：
   加载fit后，确保原始数据（data_clean）仍存在于环境中（绘图时需要分组变量IL5_group）。
   如果数据可能变动，建议将数据一并保存：
   save(fit, data_clean, file = "survival_analysis_objects.RData")
   跨平台兼容性：
   .RData和.rds文件是跨平台的，可在Windows/macOS/Linux间共享。
   版本兼容性：
   如果R版本差异较大，加载时可能出现警告，建议在相同版本环境中使用。

保存绘制的生存曲线（如图片或PDF）

1. 使用 survminer (ggplot2-based) 保存图片
   如果你用ggsurvplot()绘制生存曲线（基于ggplot2），可以用ggsave()或print()保存为图片或PDF。

方法1：直接保存为图片/PDF

library(survminer)
library(ggplot2)# 绘制生存曲线
p <- ggsurvplot(fit, data = data_clean,pval = TRUE, risk.table = TRUE,palette = "jco")# 保存为PNG（高分辨率，600 dpi）
ggsave("survival_curve.png", plot = p$plot,width = 8, height = 6, dpi = 600)# 保存为PDF（矢量图，适合论文）
ggsave("survival_curve.pdf", plot = p$plot,width = 8, height = 6)

补充：保存PDF时不支持中文解决方法：

install.packages("showtext")
library(showtext)
font_add("SimSun", "simsun.ttc")  # Windows 系统自带宋体
showtext_auto()  # 自动启用 showtext
ggsave("survival_curve.pdf", plot = surv_plot$plot,width = 8, height = 6, device = cairo_pdf)

方法2：保存带风险表的完整图形
ggsurvplot()返回的p是一个列表，包含生存曲线（p$plot）和风险表（p$table）。若需保存两者组合：

# 将图形和风险表组合
combined_plot <- arrange_ggsurvplots(list(p), print = FALSE)# 保存组合图
ggsave("survival_with_risk_table.png", combined_plot,width = 10, height = 8, dpi = 300)

2. 使用基础R图形保存
   如果使用基础R的plot()函数，可以用以下方式保存：

保存为PNG/PDF

# 打开图形设备（PNG）
png("survival_curve_baseR.png", width = 800, height = 600, res = 150)
plot(fit, col = c("red", "blue"), xlab = "Time (months)", ylab = "Survival Probability")
legend("topright", legend = levels(data_clean$IL5_group),col = c("red", "blue"), lty = 1)
dev.off()  # 关闭设备并保存# 保存为PDF（矢量图）
pdf("survival_curve_baseR.pdf", width = 8, height = 6)
plot(fit, col = c("red", "blue"), xlab = "Time (months)", ylab = "Survival Probability")
legend("topright", legend = levels(data_clean$IL5_group),col = c("red", "blue"), lty = 1)
dev.off()

3. 保存为可编辑的矢量图（如EPS/SVG）
   适合进一步在Illustrator或Inkscape中编辑：

通过ggsave保存为SVG

library(svglite)
ggsave("survival_curve.svg", plot = p$plot, device = svglite,width = 8, height = 6)

通过基础R保存为EPS

setEPS()
postscript("survival_curve.eps", width = 8, height = 6)
plot(fit, col = c("red", "blue"), xlab = "Time (months)", ylab = "Survival Probability")
legend("topright", legend = levels(data_clean$IL5_group),col = c("red", "blue"), lty = 1)
dev.off()

4. 注意事项
   文件路径：
   如果不指定路径，文件会保存在当前工作目录（通过getwd()查看）。
   自定义路径示例：
   ggsave("C:/Users/Name/Desktop/survival_curve.png", plot = p$plot)
   图形分辨率：
   调整dpi（每英寸点数）和width/height控制清晰度，尤其是PNG格式。
   多图形保存：
   如果需要保存多个分面图形（如不同亚组），用循环或lapply：
   plots_list <- lapply(subgroups, function(group) {
   ggsurvplot(subset_fit, data = subset_data, …)
   })
   for (i in seq_along(plots_list)) {
   ggsave(paste0(“survival_group_”, i, “.png”), plot = plots_list[[i]])