审稿人质疑生存分析分组不正确？这是我见到最有效的解决办法！（附代码）

通过best Separation快速获取最佳生存分组

大家好，我是阿琛。在生存分析过程中，如何对患者进行分组，以获得一个最佳的分析结果，往往是一个值得考虑的问题。在常见的分析过程中，一般存在两种选择：第一种是根据表达的中位值进行分成高表达组和低表达组；第二种则是使用最佳的截断值，即cut-off值，获得最低的P值。

然而，尽管第一种分组方法相对比较简单，但是其最终的分析结果往往并不理想，无法得到有意义的结果。因此，根据最佳的截断值来获得最佳的结果是十分必要的。

下面，我们来看下两种分组方法的分析过程。

1. R包和数据准备

1.1 R包安装与读取

#镜像设置options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")install.packages("survival")install.packages("survminer")在此，我们分别使用survival包和survminer包来构建模型和可视化展示。library(survival)library(survminer)安装完成后，使用library()函数读取相应的R包。

1.2 数据读取

rt <-read.table("gene_Surv.txt", sep="\t", header = T,check.names = F,row.names = 1)str(rt)

对于输入数据，在整理过程中，至少包含三个信息，分别是基因表达水平、患者的生存时间、生存状态等。

cli <- rt[order(rt$expression),]head(rt)

随后，使用order()函数对数据框的行按照基因表达水平从低到高进行排序。

2. 根据中位数进行分组分析

首先，我们来回顾一下如何根据表达的中位值进行分组，然后绘制相应的生存曲线。

2.1 计算P值

ssdf <- cbind(cli, data.frame(type = ifelse(cli$expression > median(cli$expression),"High","Low")))

在正式分析前，根据表达的中位值，分别将患者分成高表达组（High）和低表达组（Low）。

ssdf$type <- factor(ssdf$type, levels = c("Low","High")) # 修改因子等级table(ssdf$type)

同时，将type类型设置为因子形式，并设置相应的因子等级。注意，将分组转换为因子，以防在后续绘制生存曲线时将两条曲线弄反。使用table()函数可以看到，其中高表达患者和低表达患者各135例。

接下来，我们来进行生存分析。首先，我们需要使用 survival 包的survdiff()函数来进行统计计算。

fitd <-survdiff(Surv(futime,fustat) ~ type,data = ssdf,na.action = na.exclude)fitd

结果显示，经过生存分析，根据中位值分组后，通过 fitd 可以查得最终的P值为0.1。

pValue <- 1 - pchisq(fitd$chisq, length(fitd$n) - 1)pValuep.lab <- paste0("P", ifelse(pValue < 0.001, " < 0.001", paste0(" = ",round(pValue, 3))))

进一步计算确切的P值；结果显示，P值等于0.112。

2.2 拟合生存曲线

fit <- survfit(Surv(futime,fustat) ~ type, data = ssdf)

接下来，使用survival包的survfit()函数拟合生存曲线。survfit()函数根据先前拟合的模型来创建生存曲线，以用于后续生存曲线绘制时的输入文件。

2.3 绘制生存曲线

ggsurvplot(fit, data = ssdf, conf.int = F, pval = p.lab, surv.median.line = "hv")

使用survminer包的ggsurvplot()函数绘制生存曲线，其中将上一步拟合得到的fit模型作为输入，使用数据集ssdf进行分析。同时，设置参数conf.int = F表示图中不添加可信区间，参数surv.median.line = "hv"表示添加中位生存时间线。

3. 根据最佳cut-off值进行分组分析

计算完中位值分组的结果后，接下来我们看下如何寻找最佳的截断值。

3.1 计算不同分组截断值分组下P值

pvals<-c()hrs<-c()for(n in 1:(nrow(cli)-1)){ ssdf <- cbind(cli, data.frame(type = rep(c("Low","High"), c(n, nrow(cli)-n)))) ssdf$type <- factor(ssdf$type, levels = c("Low","High")) fitd <- survdiff(Surv(futime,fustat) ~ type, data=ssdf, na.action = na.exclude) pValue <- 1 - pchisq(fitd$chisq, length(fitd$n) - 1) pvals <- c(pvals, pValue) hr <- fitd$obs[1] * fitd$exp[2]/(fitd$obs[2] * fitd$exp[1]) hrs <- c(hrs, hr)}

由于数据集中基因表达水平是由低向高排列，因此，通过设置循环，从1个低表达患者开始，逐个计算其对应的P值和HR值，并将计算结果进行保存。最终得到269个对应的P值和HR值。

3.2 可视化展示P值

根据分组情况，结合对应的HR值和P值，构建结果对应的数据框。

dat <-data.frame(Num = 1:(nrow(cli)-1),HR = hrs,Pvalue = pvals)head(dat)

结果显示：在这个新构建的数据集中，分别展示了在不同分组情况下对应的HR和P值；例如，当1个低表达和269个高表达的情况下，其对应的HR值为3.434，P值为0.1910。

#install.packages("ggplot2")library(ggplot2)ggplot(dat, aes(x = log2(HR), y = -log10(Pvalue))) + geom_point(shape = 21, aes(fill = Num)) + scale_fill_continuous() + geom_vline(xintercept = c(-1, 1),linetype = "dashed") + geom_hline(yintercept = (-log10(0.05)),linetype = "dashed") + theme_bw()

为了进一步展示分析结果，使用ggplot()函数进行可视化展示。

结果显示，其横坐标为log2(HR)，纵坐标为-log10(Pvalue)，而且在右上角那点时P值最小。

bestNum <- dat$Num[which.min(pvals)]bestNum

获取最小P值对应的分组情况，结果显示，当Num为6时，此时对应的P值最小，即低表达患者分组6位，高表达患者分组264位。

3.3 绘制最佳分组的生存曲线

接着，根据分析得到的最佳分组，来进行生存曲线的绘制。

ssdf2 <- cbind(cli, data.frame(type = rep(c("Low","High"), c(bestNum, nrow(cli)-bestNum))))ssdf2$type <- factor(ssdf2$type, levels = c("Low","High"))table(ssdf2$type)

结果显示，此时低表达患者6位，高表达患者264位。

fitd2 <- survdiff(Surv(futime,fustat) ~ type, data = ssdf2, na.action = na.exclude)fitd2pValue2 <- 1 - pchisq(fitd2$chisq, length(fitd2$n) - 1)pValue2 #查看P值p.lab2 <- paste0("P", ifelse(pValue2 < 0.001, " < 0.001", paste0(" = ",round(pValue, 3))))####拟合生存曲线fit2 <- survfit(Surv(futime,fustat) ~ type, data = ssdf2)

同样，构建生存模型和拟合生存曲线，结果显示P值为0.02。