TCGA数据下载以及生存分析

0、R包准备
先加载所需要的R包（提前下载好）

# 加载R包
library(TCGAbiolinks)
library(tibble)
library(SummarizedExperiment)
library("stringr")
library(dplyr)
library(stats4)
library(BiocGenerics)
library(parallel)
library("AnnotationDbi")
library("org.Hs.eg.db")
library(survival)
library(survminer)

检查下TCGAbiolinks版本(如果版本过低可能会导致TCGA数据读取失败)

packageVersion("TCGAbiolinks")
‘2.29.6’

1、数据下载
以TCGA中肝癌的RNA数据为示例(可直接将数据下载到本地)

TCGAbiolinks:::getGDCprojects()$project_id
cancer_type="TCGA-LIHC"
clinical=GDCquery_clinic(project=cancer_type,type="clinical")
dim(clinical)
clinical[1:4,1:4]
head(colnames(clinical))
query <- GDCquery(project =cancer_type,
                   data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification",
                  workflow.type = "STAR - Counts"
                    )

GDCdownload(query, method = "api", files.per.chunk = 50)
expdat <- GDCprepare(query)###读取下载好的数据

2、数据整理

counts <- as.data.frame(assay(expdat))#默认提取counts数据
head(counts)###矩阵预览

counts数据

group_list=ifelse(as.numeric(substr(colnames(counts),14,15)) < 10,'tumor','normal')
table(group_list)
#仅取tumor样本
exp_tumor=counts[,group_list=='tumor']

为什么小于10的是tumor，大于10的是normal呢，因为Sample号从01-29的，其中01-09是tumor，也就是癌症样本；其中10-29是normal，也就是癌旁

TCGA编号

meta = clinical
#第一列作为行名
meta <- column_to_rownames(meta,var = "submitter_id")
colnames(meta)
#筛选我们感兴趣的临床信息
meta=meta[,colnames(meta) %in% c("vital_status",
                                 "days_to_last_follow_up",
                                 "days_to_death",
                                 "race",
                                 "gender",
                                 "age_at_index",
                                 "tumor_stage")]

#调整、筛选临床样本信息数量与顺序与表达矩阵相同
meta=meta[match(substr(colnames(exp_tumor),1,12),rownames(meta)),]

计算生存时间

#
#days_to_last_follow_up：距离死亡的时间
#day_to_last_follow_up：最后一个随访时间到开始时间
meta$days_to_death[is.na(meta$days_to_death)] <- 0   #缺失值标记为0
meta$days_to_last_follow_up[is.na(meta$days_to_last_follow_up)] <- 0
head(meta)

mata数据预览

#时间以月份记，保留两位小数
meta$time=round(meta$days/30,2)
meta$time = as.numeric(meta$time)
#2、根据生死定义活着是0，死的是1
meta$event=ifelse(meta$vital_status=='Alive',0,1)
table(meta$event)

#3 年龄分组(部分样本缺失，考虑可能的影响应该不大)
meta$age_at_index[is.na(meta$age_at_index)] <- 0
meta$age_at_index=as.numeric(meta$age_at_index)
meta$age_group=ifelse(meta$age_at_index>median(meta$age_at_index),'older','younger')

3、将ENS号转换为基因名

#####将ENS号转换为基因名
ENSEMBL=unlist(str_split(rownames(exprSet),"[.]",simplify=T))[,1]####去除ENS的版本号
df<-as.data.frame(ENSEMBL)
df$symbol <- mapIds(org.Hs.eg.db,

                     keys=ENSEMBL,

                     column="SYMBOL",

                     keytype="ENSEMBL")

转换完成

将基因名添加到exprSet数据中并作为行名

exprSet$symbol<-df$symbol
exprSet <- exprSet[!duplicated(exprSet$symbol)& !is.na(exprSet$symbol), ]####去除重复并且删掉NA
rownames(exprSet)<-exprSet$symbol
exprSet <- dplyr::select(exprSet,-symbol)

将基因名作为数据行名

接下来选择目标基因进行生存分析即可

exprSet2<-as.matrix(exprSet)###不转矩阵无法进行后续的计算
#######选择目标基因进行生存分析 
meta$gene = ifelse(exprSet2["MS4A1",]>median(exprSet2["MS4A1",]),'high','low')
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)

生存分析