R 语言小白绘图系列
第 36 弹|双基因生存曲线
1 双基因生存曲线简介
双基因生存曲线是一种用于分析生存数据的统计方法，通常在生物医学和生物统计学中使用。以下是关于双基因生存曲线的简要介绍：
定义
双基因生存曲线用于评估两个基因（或变量）对生存时间的联合影响，通常以Kaplan-Meier曲线的形式呈现。
主要步骤
数据收集：
收集相关生存数据，包括生存时间、状态（生存或死亡）、以及两个基因的表达水平。
分组：
根据基因表达水平将样本分为不同组（如高表达组和低表达组）。
生存分析：
使用Kaplan-Meier方法计算各组的生存曲线，并进行对比。
统计检验：
进行Log-rank检验等统计检验，评估不同组之间的生存差异是否显著。
应用
临床研究：分析特定基因对疾病预后的影响。 个体化医疗：根据基因表达情况，制定更合适的治疗方案。
这种方法的具体实施通常需要使用统计软件（如R或Python）来处理数据并生成生存曲线。
应用实例
如下图，根据两个基因的高低表达，将样本分为四组（高-高，低-低，高-低，低-高），然后比较组间生存是否具有差异。
横坐标代表生存时间，单位是年；纵坐标代表生存率；随着时间推移，可以看到生存率逐渐下降。通过P值可以了解组间生存率是否存在差异，若p<0.05则说明组间的生存率存在差异；
如下图，p<0.001,说明组间生存具有差异。同时可以看到，两基因均高表达时，生存是最差的；两基因均低表达时，生存是最优的。
2 源文件
输入文件有五列信息
id: 样品名称 futime:生存时间，单位为年。 fustat:生存状态，0代表存活，1代表死亡。 VCAN基因表达量 CXCR4基因表达量
3 代码
环境准备
这段代码的目的是进行包含两个变量的生存分析，以下是逐行解释：
# 安装和加载需要的R包
# install.packages("survival")
# install.packages("survminer")
library(survival) # 加载生存分析包
library(survminer) # 加载生存图形化包
这部分注释掉的代码是用于安装survival和survminer两个包，后面两行代码是加载这两个包，以便使用其中的函数进行生存分析和可视化。
inputFile = "input.txt" # 输入文件
outFile = "survival.pdf" # 输出PDF文件
var1 = "VCAN" # 用于生存分析的变量1
var2 = "CXCR4" # 用于生存分析的变量2
setwd("D:\\biowolf\\bioR\\37.survival2vars") # 设置工作目录
这部分代码定义了几个重要的变量：
inputFile: 指定输入数据文件的名称为"input.txt"。outFile: 指定输出文件的名称为"survival.pdf"。var1和var2: 这两个变量分别代表用于生存分析的生物标志物（或基因），在后续分析中将被使用。setwd(): 设置当前的工作目录，以便后续读取输入文件和保存输出文件。
rt = read.table(inputFile, header=T, sep="\t", check.names=F) # 读取输入数据
这一行代码使用read.table函数读取指定的输入文件，参数含义如下：
header=T: 表示文件的第一行包含列名。sep="\t": 指定文件使用制表符（Tab）分隔。check.names=F: 表示在读取列名时不进行修改，保留原有的列名。
总结
整段代码的功能是准备进行包含两个变量（VCAN和CXCR4）的生存分析，设置输入输出文件和工作目录，并读取输入数据。后续分析会基于这些准备工作进行。
根据基因表达，对数据分组
这段代码的作用是根据两个基因（var1和var2）的表达水平将数据分组，具体解释如下：
# 根据基因表达，对数据分组
a = ifelse(rt[, var1] < median(rt[, var1]), paste0(var1, " low"), paste0(var1, " high"))
b = ifelse(rt[, var2] < median(rt[, var2]), paste0(var2, " low"), paste0(var2, " high"))
Type = paste(a, "+", b)
rt = cbind(rt, Type)
分组操作：
类似地，对 var2的表达值进行比较，结果存储在b中。这行代码对 var1的表达值进行比较，若小于中位数则标记为"VCAN low"，否则标记为"VCAN high"。结果存储在a中。ifelse(rt[, var1] < median(rt[, var1]), paste0(var1, " low"), paste0(var1, " high"))：ifelse(rt[, var2] < median(rt[, var2]), paste0(var2, " low"), paste0(var2, " high"))：
组合分组信息：
将 a和b的分组信息结合，形成一个新的分组标签，如"VCAN low + CXCR4 high"。Type = paste(a, "+", b)：
将分组信息添加到数据框中：
将新的 Type列添加到原始数据框rt中，形成一个新的数据框，包含基因表达分组的信息。rt = cbind(rt, Type)：
总结
这段代码通过比较基因表达值与中位数，创建了一个新的分组变量Type，用于后续的生存分析，便于分析不同组合的生存情况。
生存差异统计
这段代码的目的是进行生存差异统计，具体步骤如下：
# 生存差异统计
length = length(levels(factor(Type))) # 计算分组的数量
diff = survdiff(Surv(futime, fustat) ~ Type, data = rt) # 进行生存差异分析
pValue = 1 - pchisq(diff$chisq, df = length - 1) # 计算p值
if (pValue < 0.001) {
pValue = "p<0.001" # 如果p值小于0.001，则设定为"p<0.001"
} else {
pValue = paste0("p=", sprintf("%.03f", pValue)) # 否则格式化p值
}
fit <- survfit(Surv(futime, fustat) ~ Type, data = rt) # 拟合生存曲线
逐行解释：
计算分组数量：
使用 factor()函数将Type转换为因子类型，levels()获取因子的水平（不同分组），然后用length()计算这些分组的数量。length = length(levels(factor(Type)))：
进行生存差异分析：
使用 survdiff()函数进行生存分析，比较不同Type分组的生存差异，返回的结果包含了卡方统计量等信息。diff = survdiff(Surv(futime, fustat) ~ Type, data = rt)：
计算p值：
使用卡方分布的累积分布函数 pchisq()计算p值，diff$chisq是前一步得到的卡方统计量，df为自由度（分组数减1）。pValue = 1 - pchisq(diff$chisq, df = length - 1)：
格式化p值：
if语句判断p值是否小于0.001，如果是，则赋值为"p<0.001"；否则格式化p值为三位小数的字符串。
拟合生存曲线：
使用 survfit()函数根据生存时间futime和生存状态fustat拟合生存曲线，按Type分组进行。fit <- survfit(Surv(futime, fustat) ~ Type, data = rt)：
总结
这段代码的核心是通过生存分析方法比较不同分组的生存差异，并计算相应的p值，最终得到生存曲线的拟合结果。
绘制生存曲线
这段代码的目的是绘制生存曲线并保存为PDF文件，具体步骤如下：
# 绘制生存曲线
surPlot = ggsurvplot(fit,
data = rt,
conf.int = F, # 不包含置信区间
pval = pValue,
pval.size = 5,
legend.labs = levels(factor(rt[,"Type"])), # 图例标签
legend.title = "Type", # 图例标题
legend = "right", # 图例位置
xlab = "Time(years)", # x轴标签
break.time.by = 1, # x轴刻度间隔
risk.table.title = "", # 风险表格标题
risk.table = F, # 不显示风险表
risk.table.height = .25) # 风险表高度
pdf(file = outFile, onefile = FALSE, width = 7.5, height = 5) # 设置PDF输出
print(surPlot) # 打印生存曲线
dev.off() # 关闭PDF设备
逐行解释：
绘制生存曲线：
使用 ggsurvplot()函数绘制生存曲线，传入拟合结果fit和数据集rt。surPlot = ggsurvplot(...)：
参数设置：
conf.int = F：不显示置信区间。pval = pValue：显示生存差异的p值。pval.size = 5：设置p值的字体大小。legend.labs = levels(factor(rt[,"Type"]))：根据Type的水平设置图例标签。legend.title = "Type"：图例的标题。legend = "right"：设置图例在图的右侧。xlab = "Time(years)"：x轴的标签为“Time(years)”。break.time.by = 1：x轴刻度每1单位。risk.table.title = ""：设置风险表格的标题为空。risk.table = F：不显示风险表格。risk.table.height = .25：如果显示风险表格，则设置高度为0.25。
保存为PDF：
开始一个新的PDF设备，将生存曲线输出到指定的文件，设置宽度为7.5，高度为5。 pdf(file = outFile, onefile = FALSE, width = 7.5, height = 5)：
打印生存曲线：
print(surPlot)：将绘制的生存曲线打印到PDF文件中。
关闭PDF设备：
dev.off()：关闭当前的PDF设备，完成文件的写入。
总结
这段代码的核心功能是绘制生存曲线，并配置相应的图例、标签和输出选项，最终将图形保存为PDF文件。
