生存曲线一个循环可以画太太太太太太多多多多多多图了吧

一文学会基因突变生存曲线的批量绘制

各位解螺旋的小伙伴们大家周日好,我是先锋宇,欢迎大家来到每周日的先锋宇专栏。上一期的推文我们打开了基因组学的大门,不知道大家对于好看又好吃的棒棒糖图有没有留下深刻的印象呢?棒棒糖虽好,但是也不要贪吃哦。今天小编继续给大家带来另外一份大餐,我们来进行生存曲线再一次的批量绘制,这一次我们把目光从转录组学移到基因组学,我们在基因组学进行基因突变生存曲线的批量绘制,小伙伴们有没有很心动呢?那怀着愉悦的心情和小编一起来学习一下吧

写在前面

生存曲线的绘制对于生信分析是非常重要的,往往在一篇文章中我们看到肯定不止绘制一条生存曲线,而是有很多条生存曲线。先锋宇老师很喜欢说一句话:机械重复的工作拿给机器来完成。所以我们应该摆脱机械重复的工作,而应该用更多的时间来享受生活。那么今天我们就一起来看看基因组学的生存曲线批量绘制能不能提升我们的科研效率呢?

代码演示

经过上一期推文的讲解,我们已经从Xena下载了胃癌的基因组学数据以及过滤掉了不引起氨基酸改变的突变,得到了如下的数据框(这一步不清楚的同学可以去复习一下先锋宇老师上一期的推文:小白重头学习基因组学分析,从绘制棒棒糖图开始)


那接下来为了绘制每个基因突变与否的批量生存曲线,我们首先需要去构建突变矩阵,指明基因在哪些样本是突变的,哪些样本是不突变的。所以我们首先需要提取上面这个数据框的前两列(这里不懂突变矩阵没关系,看完了你就都明白了)


mutation_matrix <- verscan2_code %>% dplyr::select(1,2)

可以看到这两列数据是一个明显的长数据,里面只有样本对应有突变的基因,所以我们先添加一列mutation,在后面标记上这些都是样本里面的有突变的基因。


mutation_matrix$group <- c("mutation")

可以看到我们添加了一列信息来指明这些样本中的基因是突变的,这在我们接下来构建突变矩阵的时候非常有用。接下来我们需要对突变矩阵去重,因为一个样本的同一个基因发生不同的突变我们在突变矩阵里面就只需要知道它是mutation,不需要管具体的突变方式。所以我们这里使用unique函数进行去重复。


mutation_matrix <- unique(mutation_matrix)

可以看到去重之后我们去掉了20000多行,那么接下来我们就需要把这个长数据变成宽数据,也就是我们习惯的矩阵的形式,这样才便于我们进行作图。长宽数据转换在R里面用到的也比较多,还不会的小伙伴们也可以参考我之前的推文,我详细介绍了长宽数据转换的三种方法,这里我们使用比较常用的pivot_wider函数进行转换,函数作用和名字一样,把数据变wider。


mutation_matrix1 <- pivot_wider(mutation_matrix,                               names_from = "Sample_ID",                               values_from = "group")

可以看到变为宽数据之后,列名为样本名,然后我刚刚加的一列标签(mutation)在这里就起作用啦,当数据变为宽数据之后,我们就能看到哪些基因在哪里样本里面是突变的,而没有突变的在这里就是NA,那么接下来我们再把NA全部变为wild(野生型)就可以啦。


mutation_matrix1 <- mutation_matrix1 %>% column_to_rownames("gene")mutation_matrix1[is.na(mutation_matrix1)] <- c("wild")

可以看到突变矩阵就构建好了,那么接下来我们为了要绘制生存曲线,我们接下来就需要把TCGA样本的生存数据引入进来,前面先锋宇专栏的推文已经讲解了生存数据的下载和读取,这里就不加赘述啦,我们接下来读取TCGA胃癌的生存数据。


stad_survival <- read_tsv("TCGA-STAD.survival.tsv.gz") %>% dplyr::select(1,2,4)colnames(stad_survival) <- c("sample", "fustat", "futime")stad_survival$futime <- stad_survival$futime/365

可以看到我们把TCGA胃癌病人的生存信息读入了,那么接下来我们就把突变矩阵和生存信息进行整合,根据TCGA样本讲两个数据框整理为一个数据框,我们这里仍然不使用merge函数,我们使用更为简单的inner_join函数来进行操作


mutation_matrix2 <- t(mutation_matrix1) %>% as.data.frame() %>% rownames_to_column("sample")mutation_survival <- inner_join(stad_survival, mutation_matrix2,                               by = "sample")

可以看到我们把胃癌病人的生存信息与基因突变信息整合起来,这样接下来我们就可以进行基因突变生存曲线的批量绘制啦。

首先为了降低难度,我们仍然先进行一条生存曲线的绘制,所以我们提取前面4列,只保留矩阵的第一个基因和生存信息。


gene_survival <- mutation_survival %>% dplyr::select(1,2,3,4)

接下来我们进行生存曲线的绘制,我们首先打开survminer和survival包,方便我们调用survdiff函数来计算pvalue。


library(survminer)library(survival)diff=survdiff(Surv(futime, fustat) ~ gene_survival[[4]],data = gene_survival)pValue=1-pchisq(diff$chisq,df=1)

接下来我们使用survfit来构建一个生存函数,并使用ggsurvplot函数来进行生存曲线绘制。


fit <- survfit(Surv(futime, fustat) ~ gene_survival[[4]],data = gene_survival)ggsurvplot(fit,           data=gene_survival,          conf.int=TRUE,          pval=pValue,          pval.size=6,          legend.labs=c("mutation", "wild"),          legend.title="CLSTN1",          xlab="Time (years)",          ylab="Overall survival",          break.time.by = 1,          risk.table.title="",          palette=c("#d7191c", "#1a9641"),          risk.table=T,          risk.table.height=.25)

可以看到一张可以用于发表的基因组的生存曲线就绘制好啦。

进阶

一张生存曲线绘制好了,但是这还远远不敢,成年人的世界里我们不做选择,而是全部都要,所以我们就要对生存曲线进行批量出图,其实总的思想跟前面讲的转录组也是一样的,主要的区别就是这里不是基因的表达量,而是基因的突变信息,同样是按照列对每个基因分别进行生存曲线的绘制。


for (gene in colnames(mutation_survival[4:17413])) { diff=survdiff(Surv(futime, fustat) ~ mutation_survival[[gene]],data = mutation_survival) pValue=1-pchisq(diff$chisq,df=1) fit <- survfit(Surv(futime, fustat) ~ mutation_survival[[gene]],data = mutation_survival) surPlot = ggsurvplot(fit,             data=gene_survival,            conf.int=TRUE,            pval=pValue,            pval.size=6,            legend.labs=c("mutation", "wild"),            legend.title=gene,            xlab="Time (years)",            ylab="Overall survival",            break.time.by = 1,            risk.table.title="",            palette=c("#d7191c", "#1a9641"),            risk.table=T,            risk.table.height=.25) pdf(file=paste0("survival/",gene,".pdf"), onefile=FALSE, width=6.5, height=5.5) print(surPlot) dev.off()}

这里我们使用for循环对于列进行循环,最终我们可以得到10000多个基因突变信息的生存曲线图,大家只需要慢慢挑选自己想要的基因的生存曲线图进行展示即可。


写在最后

机械重复的工作我们应该尽可能交给计算机去完成,如果我们每天就是去做机械重复的工作,比如这里我们去重复几百次甚至上万次一样的操作去绘制想要的生存曲线,那这样科研的效率就会大打折扣,有些时候学会适当的偷懒,把这些时间用来享受生活,这才是我们应该经历的人生,而不是永远沉浸在科研的无底洞里面,看不到尽头。


好啦,今天的生存曲线批量绘制就到这里结束啦,希望这道大餐能给冬日冷冷的天气里面增加一丝暖意。我是先锋宇,我们下周日继续相约挑圈联靠的先锋宇专栏,我们不见不散。

往期传送门:

让生信工作者失业的神器——DrBioRight,真舍不得告诉你!

高效数据清洗!这个R包太强大了!你一定要试试!(文末附赠小彩蛋)

一站式分析R包来了,承包了生信各种分析!太全能了!

搞定这一步,说明你学R有天赋!TCGA数据从下载到差异分析(附代码)

别走,baby,COX森林需要你

生存曲线居然能够批量绘制了

ROC居然能够批量筛选

TCGA转录组和临床信息联合分析,结果发现......

小白重头学习基因组学分析,从绘制棒棒糖图开始

欢迎大家关注解螺旋生信频道-挑圈联靠公号~


END
撰文丨先锋宇
排版丨四金兄
值班 | 菠小萝
主编丨小雪球



2021年国自然申请指南已公布

2021申请国自然还很迷茫?
国自然申请指南找不到重点?

酸谈为你解读最新国自然申请指南!

本周四谈老谈直播主题
---2021年国自然申请指南解读
---

下方视频号点击预约直播
  酸菜老助你冲刺国自然!

👇👇👇

周四酸谈直播基本信息

主题:《2021年国自然申请指南解读》
日期:1月21日 18:00—20:00
(👆点击“预约”就完事了!!!)

生物医学科研方法

癌症新发现:这种肠道细菌能够转移到乳腺中,诱导乳腺增生、促进乳腺癌

2021-1-18 6:14:47

生物医学科研方法

河北高校放假后共有2127名大学生滞留石家庄

2021-1-18 6:34:44