如何保证自己的生存分析结果图有意义
介绍
一般来说，我们做生存分析，会有（P<0.05）和（P>0.05）两种结果。KM plot在生物医学中很常见，主要用来做预后分析，比如可以根据表达量把病人分成两组，然后比较哪组病人预后好，进而可以得出基因表达量高低与病人预后好坏相关性的结论。 画KM plot时，有时候会比较纠结怎样对病人进行分组，如何来设置分组的cutoff。一般来说常见的几种设置cutoff值得思路如下： 1：大多数情况下，根据表达量从低到高对样本进行排序，取前50%为低表达，后50%为高表达，然后画KM plot。 2：还有一些文章也会将样本表达量均分为三组或者四组。 3：一些文章也会选一些其它的cutoff，比如前1/3和后2/3，前25%和后25%（中间50%的数据去掉）。
例子
例如下面例子所示：（通过NFE2L2基因的表达量中位值，我们将所有的样本分为高表达和低表达两组，然后通过绘制KM生存分析曲线的形式来探讨两组生存概率是否存在差别）
> # =======================================================
> ####
> # =======================================================
>
> library(dplyr)
>
> library(tidyr)
>
> library(tidyverse)
>
> library(survival)
>
> library(survminer)
>
> setwd('D:\\train\\sur')
>
> rm(list=ls())
>
> data <- read.csv('data.csv', header = T)
>
> head(data)
sample NFE2L2 OS OS.Time
1 sam590 660 0 1.413699
2 sam1034 749 0 2.326027
3 sam164 859 0 1.227397
4 sam1079 1034 0 5.569863
5 sam266 1068 0 1.221918
6 sam935 1098 0 1.986301
>
> str(data)
'data.frame': 1086 obs. of 4 variables:
$ sample : Factor w/ 1086 levels "sam1","sam10",..: 633 41 160 90 273 1016 32 37 271 400 ...
$ NFE2L2 : int 660 749 859 1034 1068 1098 1126 1197 1201 1208 ...
$ OS : int 0 0 0 0 0 0 0 0 0 0 ...
$ OS.Time: num 1.41 2.33 1.23 5.57 1.22 ...
>
> data['NFE2L2'] <- ifelse(data['NFE2L2'] > median(data[,'NFE2L2']), 'high', 'low' )
>
> fit <- survfit(Surv(OS.Time, OS) ~NFE2L2, data =data)
>
> ggpar( ggsurvplot(
+ fit, data = data,
+ pval = TRUE,
+ # conf.int = TRUE,
+ # legend.title = "CUL3",
+ xlim = c(0,5),
+ xlab = "Time in years",
+ break.time.by = 1,
+ pval.size = 8,
+ risk.table = "absolute",
+ risk.table.y.text.col = T,
+ risk.table.y.text = FALSE,
+ risk.table.fontsize = 5,
+ # legend.labs = c("high", "Low"),
+ palette = c("#E41A1C","#377EB8")),
+ font.y = c(16, "bold"),
+ font.x = c(16, "bold"),
+ legend = "top",
+ font.legend = c(16, "bold"))
但是很尴尬的发现，通过该基因的中位值分成的高低两组并不存在生存差异。这个时候，我们首先可以通过三分法或者四分法，将患者均分为三个组别或者四个组别。
# =======================================================
####
# =======================================================
library(dplyr)
library(tidyr)
library(tidyverse)
library(survival)
library(survminer)
setwd('D:\\train\\sur')
rm(list=ls())
data <- read.csv('data.csv', header = T)
head(data)
str(data)
rt <- data %>%
mutate_at(vars(NFE2L2:NFE2L2),
funs(Hmisc::cut2(as.numeric(.), g = 3)))
table(rt$NFE2L2)
fit <- survfit(Surv(OS.Time, OS) ~NFE2L2, data =rt)
ggpar( ggsurvplot(
fit, data = rt,
pval = TRUE,
# conf.int = TRUE,
# legend.title = "CUL3",
xlim = c(0,5),
xlab = "Time in years",
break.time.by = 1,
pval.size = 8,
risk.table = "absolute",
risk.table.y.text.col = T,
risk.table.y.text = FALSE,
risk.table.fontsize = 5
# legend.labs = c("high", "Low"),
# palette = c("#E41A1C","#377EB8")
),
font.y = c(16, "bold"),
font.x = c(16, "bold"),
legend = "top",
font.legend = c(16, "bold"))
我们通过Hmisc::cut2(as.numeric(.), g = 3))，将患者均分为三组，但是虽然P值下降了，仍然没有达到想要的（P<0.05）的效果。所以，我们需要其他方法。
# =======================================================
####
# =======================================================
library(dplyr)
library(tidyr)
library(tidyverse)
library(survival)
library(survminer)
setwd('D:\\train\\sur')
rm(list=ls())
data <- read.csv('data.csv', header = T)
sur.cut <- surv_cutpoint(data, time= 'OS.Time',
event = 'OS' ,
variables = 'NFE2L2' )
sur.cut
sur.cat <- surv_categorize(sur.cut)
table(sur.cat$NFE2L2)
fit <- survfit(Surv(OS.Time, OS) ~NFE2L2, data =sur.cat)
ggpar( ggsurvplot(
fit, data = sur.cat,
pval = TRUE,
conf.int = TRUE,
# legend.title = "CUL3",
xlim = c(0,5),
xlab = "Time in years",
break.time.by = 1,
pval.size = 8,
risk.table = "absolute",
risk.table.y.text.col = T,
risk.table.y.text = FALSE,
risk.table.fontsize = 5,
# legend.labs = c("high", "Low"),
palette = c("#E41A1C","#377EB8")),
font.y = c(16, "bold"),
font.x = c(16, "bold"),
legend = "top",
font.legend = c(16, "bold"))
通过sur.cut我们达到了P小于0.05的目标，这一步的主要原理是，放弃以前所用的中位值来定义高低组的方法，采用不同的阈值来重新定义高低分组以达到最低的P值。我们看到经过这样的步骤，高组别的患者明显多于低组别。
> sur.cut
cutpoint statistic
NFE2L2 3705 3.039985sur.cut即是我们找到的最佳cutoff值，当我们通过这个cutoff值来定义高低分组时，P值达到最低。但是我们可以逐渐尝试该cutoff值附近的值，来找到一个合适的阈值。
本文参与 腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。
原始发表：2019.10.08 ，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录
腾讯云开发者
