生存分析中给基因表达量（连续变量）设置阈值

用户7010445

发布于 2020-04-30 18:14:36

1.5K1

发布于 2020-04-30 18:14:36

在论文 Construction of a Competitive endogenous RNA network and identification of potential regulatory axis in gastric cancer 中，他在生存分析的方法部分写到 A “survminer” package in R software to determine the best cut-off of the expression value for survival analysis。

通过help(package="survminer")查看帮助文档，找打了函数surv_cutpoint()和surv_categorize()

下面试着运行帮助文档中的例子数据集是这个样子滴：

> head(myeloma)
         molecular_group chr1q21_status treatment event  time
GSM50986      Cyclin D-1       3 copies       TT2     0 69.24
GSM50988      Cyclin D-2       2 copies       TT2     0 66.43
GSM50989           MMSET       2 copies       TT2     0 66.50
GSM50990           MMSET       3 copies       TT2     1 42.67
GSM50991             MAF           <NA>       TT2     0 65.00
GSM50992    Hyperdiploid       2 copies       TT2     0 65.20
           CCND1 CRIM1 DEPDC1    IRF4   TP53   WHSC1
GSM50986  9908.4 420.9  523.5 16156.5   10.0   261.9
GSM50988 16698.8  52.0   21.1 16946.2 1056.9   363.8
GSM50989   294.5 617.9  192.9  8903.9 1762.8 10042.9
GSM50990   241.9  11.9  184.7 11894.7  946.8  4931.0
GSM50991   472.6  38.8  212.0  7563.1  361.4   165.0
GSM50992   664.1  16.9  341.6 16023.4 2096.3   569.2

运行代码

library(survminer)
data("myeloma")
head(myeloma)
res.cut<-surv_cutpoint(myeloma,time="time",event="event",
                       variables = c("DEPDC1","WHSC1","CRIM1"))
summary(res.cut)
plot(res.cut,"DEPDC1",palette="npg")
res.cat<-surv_categorize(res.cut)
head(res.cat)

最后把数据集变成了：

> head(res.cat)
          time event DEPDC1 WHSC1 CRIM1
GSM50986 69.24     0   high   low  high
GSM50988 66.43     0    low   low   low
GSM50989 66.50     0    low  high  high
GSM50990 42.67     1    low  high   low
GSM50991 65.00     0    low   low   low
GSM50992 65.20     0   high   low   low

下面画图

library(survival)
fit<-survfit(Surv(time,event)~DEPDC1,data=res.cat)
ggsurvplot(fit,data=res.cat)

添加一些其他参数

ggsurvplot(fit,data=res.cat,conf.int = T,pval = T)

image.png

pval还可以传递给文本

ggsurvplot(fit,data=res.cat,conf.int = T,pval = "ABC\nDEF\nGHI")

更改图例的标题和内容用到legend.title和legned.labs

ggsurvplot(fit,data=res.cat,conf.int = T,pval = "ABC\nDEF\nGHI",
           legend.labs=c("A","B"),
           legend.title=c("D"))

再加一些其他的参数

ggsurvplot(fit,data=res.cat,conf.int = T,pval = "ABC\nDEF\nGHI",
           legend.labs=c("A","B"),
           legend.title=c("D"),
           risk.table = T,
           ncensor.plot=T,
           palette = c("#E7B800", "#2E9FDF"),
           ggtheme = theme_bw())