生信文章数据分析--TCGA生存分析

首先对该文中第一幅图进行分析，该图描述的是：通过单因素生存分析发现，在TCGA-GC患者样本中，Stromal scores与患者预后明显相关(P=0.012)。

168ee8ec78849af9a91ef6027d331c59.png

数据下载

获取肿瘤相关Scores值

Stromal scores, immune scores 和 estimate scores 是反映肿瘤细胞纯度、免疫浸润程度的一个指标，可从MD Anderson (https://bioinformatics.mdand erson.org/estimate/)下载获取，几乎包括了TCGA所有肿瘤样本的ESTIMATE score。

08e0e86bb2c377982fca5b915d89ef64.png
文章中研究的是GBM，即多形性胶质母细胞瘤，这个是下载的score数据

ID  Stromal_score   Immune_score    ESTIMATE_score
TCGA-02-0001-01 491.67  2852.31 3343.98
TCGA-02-0003-01 1248.81 3210.47 4459.27
TCGA-02-0004-01 2016.62 1566.65 3583.27

TCGA数据分析简述

此外，还需要下载获取对应的随访数据，TCGA数据的下载使用cgdsr包

对数据进行处理

#读入MD Anderson下载的数据
> estimation<-read.table("gbm_RNAseqV2.txt",                       
+                        header = T,stringsAsFactors = F)
下载获取临床数据
library("cgdsr")
library("tidyverse")

**获取GBM的临床数据**
mycgds <- CGDS("http://www.cbioportal.org/")
mycancerstudy <- "gbm_tcga"
mycaselist <- getCaseLists(mycgds,mycancerstudy)[1,1] 
myclinicaldata <- getClinicalData(mycgds,mycaselist) #共有54列信息
choose_columns <- c("AGE",'SEX','OS_STATUS','OS_MONTHS',
 "DFS_MONTHS","DFS_STATUS")
choose_clinicaldata <- myclinicaldata[,choose_columns] #只选择其中的6列用于后续分析


#合并
> estimation$ID=gsub("-",".",estimation$ID)

> dat$ID = rownames(dat)
> dat_estimation=merge(dat, estimation,by= "ID") 

分组

分组
按照文中给到的cutoff值，进行分组：2/3为higer组，1/3为lower 组。

5aac19cc5a01f60aaf16b7c67da64d2e.png

分组

> dat_estimation$Stromal_group=ifelse(  
+     dat_estimation$Stromal_score>
+         quantile(dat_estimation$Stromal_score,0.33),
+     "High","Low")

> head(dat_estimation, 3)
               ID AGE    SEX  OS_STATUS OS_MONTHS DFS_MONTHS            DFS_STATUS
1 TCGA.02.0001.01  44 Female 1:DECEASED     11.76       4.50 1:Recurred/Progressed
2 TCGA.02.0003.01  50   Male 1:DECEASED      4.73       1.31 1:Recurred/Progressed
3 TCGA.02.0004.01  59   Male 1:DECEASED     11.33      10.32 1:Recurred/Progressed
  Stromal_score Immune_score ESTIMATE_score Stromal_group
1        491.67      2852.31        3343.98          High
2       1248.81      3210.47        4459.27          High
3       2016.62      1566.65        3583.27          High

进行生存分析

> install.packages("survival")
# 加载R包l
library(survival)
# 拟合生存曲线
> my.surv <- Surv(dat_estimation$OS_MONTHS,dat_estimation$OS_STATUS=='1:DECEASED'
> kmfit2 <- survfit(my.surv~dat_estimation$Stromal_group) 
> ggsurvplot(kmfit2,data = dat_estimation, pval=T)

image.png

其他两组生存分析也是如此进行绘制。生存分析的结果与文章发表结果基本一致，Low scores组患者的生存状态要优于High scores组。