（生物信息学）R语言绘图初-中-高级——3-10分文章必备——生存曲线（初级）

楷然教你学生信

已于 2022-07-24 20:12:23 修改

阅读量3.3k

点赞数 10

分类专栏：生存曲线生物信息学与R语言 R语言作图文章标签： r语言开发语言

于 2022-07-04 21:35:08 首次发布

本文链接：https://blog.csdn.net/weixin_46500027/article/details/123245980

版权

生物信息学与R语言同时被 3 个专栏收录

10 篇文章 28 订阅

订阅专栏

生存曲线

2 篇文章 5 订阅

订阅专栏

R语言作图

2 篇文章 1 订阅

订阅专栏

生物信息学文章的发表要求除了思路和热点以外，图片绘制是否精美也是十分重要的，本专栏为（生物信息学）R语言绘图初-中-高级——3-10分文章必备，主要通过大量文献，总结3-10分文章中高频出现的各种图片，并给大家提供图片复现的R语言代码，及图片识读。

本专栏将向大家介绍的图片绘制如下：

1. 散点图

2. 箱线图

3.条形图

4.正负条形图

5.区组条形图

6.小提琴图

7.热图

8.Venn图

9.生存曲线

10.森林图

11.TSNE

12.瀑布图

13.ROC曲线

14.点阵图

15.相关系数图

16.饼图

17.树形图

18.气泡图

19.火山图

20.点图

上次我们绘制了Venn图：

（生物信息学）R语言绘图初-中-高级——3-10分文章必备——Venn图（韦恩图）（初级）_楷然教你学生信的博客-CSDN博客（生物信息学）R语言绘图初-中-高级——3-10分文章必备——Venn图（韦恩图）（初级）https://blog.csdn.net/weixin_46500027/article/details/125496246?spm=1001.2014.3001.5501

下面绘制生存曲线：

首先准备数据，数据格式如下：

这是来自TCGA—BLCA的数据，OS.time生存时间，OS生存状态：1表示死亡，0表示存活，NSUN6是随意筛选的某个基因的表达量。下面开始测试代码：

1.设置工作路径，读取数据

setwd("D:\\data")
dir()
data <- read.csv("survival.csv",header = T,sep = ",")
head(data)##查看数据
> head(data)
            Sample OS.time OS    NSUN6
1 TCGA.2F.A9KO.01A     734  1 4.015870
2 TCGA.2F.A9KP.01A     364  1 3.922716
3 TCGA.2F.A9KQ.01A    2886  0 4.120226
4 TCGA.2F.A9KR.01A    3183  1 2.909812
5 TCGA.2F.A9KT.01A    3314  0 3.160342
6 TCGA.2F.A9KW.01A     254  1 2.203534

2.对需要分析的数据进行分组，这里以NSUN6中位值为界限，高于中位值的设置为高表达组，低于中位值的设置为地表达组。当然，我们也可以用一些软件确定自变量的cut-off值，如x-tile。x-tile的使用教程将在下期分享。

data$a <- ifelse(data$NSUN6>median(data$NSUN6),"High NSUN6 expression","Low NSUN6 expression")
> head(data)
            Sample OS.time OS    NSUN6                     a
1 TCGA.2F.A9KO.01A     734  1 4.015870 High NSUN6 expression
2 TCGA.2F.A9KP.01A     364  1 3.922716 High NSUN6 expression
3 TCGA.2F.A9KQ.01A    2886  0 4.120226 High NSUN6 expression
4 TCGA.2F.A9KR.01A    3183  1 2.909812  Low NSUN6 expression
5 TCGA.2F.A9KT.01A    3314  0 3.160342 High NSUN6 expression
6 TCGA.2F.A9KW.01A     254  1 2.203534  Low NSUN6 expression

3.接下来绘制生存曲线，代码如下


library(coin)
library(survminer)
library(survival)
fit <- survfit(Surv(OS.time,OS) ~ data$a, data = data)
d <- ggsurvplot_list(fit,data = data,pval = T,##是否添加P值
                     conf.int = F,### 是否添加置信区间
                     risk.table = T, # 是否添加风险表
                     risk.table.col = "strata", 
                     ###linetype = "strata",
                     surv.median.line = "hv", # 是否添加中位生存线
                     risk.table.y.text.col = F,risk.table.y.text = FALSE,
                     ggtheme = theme_bw()+theme(legend.text = element_text(colour = c("red", "blue")))
                     +theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
                     +theme(plot.title = element_text(hjust = 0.5,size = 16,face = "bold"),axis.title.y.left = element_text(size = 16,face = "bold",vjust = 1),axis.title.x.bottom = element_text(size = 16,face = "bold",vjust = 0))
                     +theme(axis.text.x.bottom = element_text(size = 12,face = "bold",vjust = -0.8,colour = "black"))
                     +theme(axis.text.y.left = element_text(size = 12,face = "bold",vjust = 0.5,hjust = 0.5,angle = 90,colour = "black"))
                     +theme(legend.title = element_text(face = "bold",family = "Times",colour = "black",size = 12))
                     +theme(legend.text = element_text(face = "bold",family = "Times",colour = "black",size =12)), # Change ggplot2 theme
                     palette = c("#bc1e5d", "#0176bd"),##如果数据分为两组，就写两种颜色，三组就写三种颜色，具体的颜色搭配参考上一期发布的R语言颜色对比表
xlim=c(10,5000),##x轴的阈值，根据随访数据的天数进行更改
xlab = "Days")##随访的时间时天，就写Days，月份就写Months

d

代码中需要更改的地方已进行注释。

运行的图片如下：

当然需要对图片进行进一步加工，如果不想要中位生存线，风险表和P值，可以在代码中修改：

修改如下：

最终得到的图如下所示：

最后我们分析一下两组之间的风险差异，使用单因素cox回归：

library(coin)
coxreg <- coxph(Surv(OS.time,OS) ~ data$a, data = data)
summary(coxreg)###exp(coef)表示HR值

结果如下：

> summary(coxreg)###exp(coef)表示HR值
Call:
coxph(formula = Surv(OS.time, OS) ~ data$a, data = data)

  n= 403, number of events= 178 
   (因为不存在，4个观察量被删除了)

                             coef exp(coef) se(coef)     z Pr(>|z|)  
data$aLow NSUN6 expression 0.2549    1.2903   0.1512 1.686   0.0919 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
data$aLow NSUN6 expression      1.29      0.775    0.9594     1.735

Concordance= 0.542  (se = 0.02 )
Likelihood ratio test= 2.86  on 1 df,   p=0.09
Wald test            = 2.84  on 1 df,   p=0.09
Score (logrank) test = 2.86  on 1 df,   p=0.09

其中exp(coef)代表风险比（HR），1.29表示NSUN6低表达组较高表达组死亡风险增加1.29倍。但是logrank test显示P值没有意义。最后我们使用PS或者AI软件对图片进行修整，填上相关信息，最终图片如下所示：