专栏首页生物信息云TCGA数据库:生存分析

TCGA数据库:生存分析

本文介绍生存分析,其实,在R中,生存分析很简单,大家在网上能找到无数的文章。利用survival包就可以。就是按照下列公式就可以完成简单的生存分析。

fit <- survfit(Surv(生存时间, 生存状态) ~ 分组, data=数据框)

我们这里就结合基因的表达量,来进行分析。

首先加载我们的数据。

options(stringsAsFactors = F)
#加载表达数据
load("F:/TCGA/HTSeq-FPKM/Rdata/data/TCGA-COAD-Exp.Rdata")
#加载临床数据
load("K:/TCGA/clinicalData/tidyAllCancerData/TCGA-COAD -Clindata.Rdata")

表达数据和临床数据,我之前已经上传到网盘

之前处理后的数据进行简单的处理,其实就是去掉正常组织的样本,再把列名变成3个字段。因为原来表达矩阵中病人的barcode长,"TCGA-AA-3662-11A-01R-1723-07",而临床数据中的只有前3段。具体关于barcode这个在前面有介绍,也可以参考文章TCGAbiolinks包介绍,里面有详细介绍。

##提取肿瘤样本的表达矩阵
exp <- transomeData[["allGenesExp"]][["Exp"]]
tumExp <- exp[,transomeData[["allGenesExp"]][["TumorBarcode"]]]
barcode <- gsub("(.*?)-(.*?)-(.*?)-.*","\\1-\\2-\\3",colnames(tumExp))
colnames(tumExp) <- barcode 
##

同样我们也要处理一下临床数据,我们之前处理的临床数据是这样的:

我们这里也需要简单处理一下。

#提取临床数据
library(dplyr)
clin <- clinData
clin <- filter(clin,follow != "--")
clin <- clin[,c("sample","survivalState","follow")]
clin[1,3]
clin$follow <- as.numeric(clin$follow)/365

其实,就是只要生存时间和生存状态这2个数据,再删除缺失值。follow除了365,单位就是年啦。

然后我们将表达矩阵与临床数据融合,因为不是每个病人的数据都是一一对应的,简单说,就是病人有表达数据,但他的临床数据就不全,我们也删除了缺失值的病人的临床数据,所以我们只需要具有临床数据又有表达数据的病人的数据。也就是取一个交集。

##融合数据
interbarcode <- intersect(barcode,clin$sample)
gene <- "LLGL2"
geneExp <- tumExp[gene,interbarcode]
geneExp <- t(geneExp) %>% as.data.frame() %>% data.matrix() %>% as.data.frame() 
geneExp$sample <- rownames(geneExp)
mergdata <- merge(clin,geneExp,by = "sample")

这里的gene我只写了一个,所以融合后是下面这样的数据。你也可以把所有基因都融合进去,这里是案例,就演示了一个。

接下来就是添加一个分组信息。

mergdata$Group <- ifelse(mergdata[,gene] > median(mergdata[,gene]),"High","Low")

我们以表达值的中位数为分界线,高于中位值为高表达,低于或等于中位值为低表达。也可以用均值。

得到上面这样的数据后,我们就可以按照刚刚的公式进行生存分析了:

######################### 生存分析
library(survival)
library(survminer)
fit <- survfit(Surv(follow, survivalState) ~ Group, data=mergdata)

绘图的话就用ggsurvplot函数。

ggsurvplot(fit,data = mergdata)

好像不是很美观,我们可以调整一下参数,比如y轴下部分很空,我们可以调整一下y轴坐标。

这样看着就好很多啦。我们在进行其他参数的调整。

ggsurvplot(fit, main = "Survival curve",
           font.main = c(16, "bold", "darkblue"),
           font.tickslab = c(12, "plain", "darkgreen"),
           legend.title = gene,
           legend.labs = c("High", "Low"),
           ylim = c(0.75,1),
           # Add p-value and tervals
           pval = TRUE,
           palette = c("#FF0000", "#2200FF"),
           ggtheme = theme_bw() 
)

如果我们要一次批量分析很多基因的高低表达与生存的关系,写一个循环,批量绘图了。

尽管本文是介绍基因表达量的生存分析,但其他的也是一样,就看你怎么分组,比如我们前面介绍SNP的数据处理后,能否做某基因突变与野生型的生存分析呢?其实都是一样的道理,其他的也是一样。照葫芦画瓢而已,大家自己去试试。

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-05-30

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R语言基础绘图教程——第4章:面积图和饼图

    DoubleHelix
  • R绘图笔记 | 散点分布图与柱形分布图

    https://docs.qq.com/sheet/DV0dxREV1YkJ0ZmVj

    DoubleHelix
  • R语言基础绘图教程——第5章:直方图和柱状图

    DoubleHelix
  • 为什么你觉得Matplotlib用起来很困难?因为你还没看过这个思维导图

    Matplotlib是一个流行的Python库,可以很容易地用于创建数据可视化。然而,设置数据、参数、图形和绘图在每次执行新项目时都可能变得非常混乱和繁琐。而且...

    HuangWeiAI
  • python(实操4):录音文件的读取、

    py3study
  • 前端模拟ajax接口

    在平常开发中,了解完需求后,前端和后端会确定页面的需要的ajax接口,及接口的细节(请求与响应的格式)。然后,前后端就可以各自开工~ (注:在本文的接口均指aj...

    Joel
  • 数据分析实战:利用python对心脏病数据集进行分析

    今天在kaggle上看到一个心脏病数据(数据集下载地址和源码见文末),那么借此深入分析一下。

    Python进阶者
  • DRY原则的一个简单实践

    我们之前就发过一篇相关的文章:https://www.cnblogs.com/powertoolsteam/p/12758496.html 其中也提到了包括DR...

    葡萄城控件
  • 算法学习之动态数组类

    让我们的数据结构可以放置“任何”数据类型 不可以是基本数据类型,只能是类对象 boolean,byte , char,short,int,long,floa...

    慕白
  • Vue响应式原理

    Vue是数据驱动视图实现双向绑定的一种前端框架,采用的是非入侵性的响应式系统,不需要采用新的语法(扩展语法或者新的数据结构)实现对象(model)和视图(vie...

    伯爵

扫码关注云+社区

领取腾讯云代金券