matlab 生存分析,如何保证自己的生存分析结果图有意义

介绍

一般来说，我们做生存分析，会有(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

>

> 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'] median(data[,'NFE2L2']), 'high', 'low' )

>

> fit

>

> 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

head(data)

str(data)

rt %

mutate_at(vars(NFE2L2:NFE2L2),

funs(Hmisc::cut2(as.numeric(.), g = 3)))

table(rt$NFE2L2)

fit

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

sur.cut

event = 'OS' ,

variables = 'NFE2L2' )

sur.cut

sur.cat

table(sur.cat$NFE2L2)

fit

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.039985

sur.cut即是我们找到的最佳cutoff值，当我们通过这个cutoff值来定义高低分组时，P值达到最低。但是我们可以逐渐尝试该cutoff值附近的值，来找到一个合适的阈值。