宝藏：TCGA任意基因任意肿瘤，随意分析（附数据和代码）

当我们拿到一个基因，肯定是想知道它的各种基本信息，今天的教程可以做到

1 任意一个肿瘤在泛癌中的表达情况

2 任意一个基因在肿瘤和正常组织中的表达情况

3 任意一个基因在肿瘤不同分期，不同性别等临床特性的表达情况

4 任意一个基因在任意一个肿瘤，或肿瘤的某种特征中的生存分析

5 这个数据能做到的太多，只要充分发挥想象力，所有的数据获取方式见文末

首先我们从UCSC Xena数据框下载pancancer的标准化后的表达谱和临床资料

放到我们的工作目录上，解压

读取数据 

这一步耗时较长（如果这个过程你的电脑hold不住了，可以直接用后面整理好的数据，开始作图）

rm(list = ls())library(tidyverse)#读取数据ALLdata <- data.table::fread('tcga_RSEM_gene_tpm',data.table = F)ALLdata[1:5,1:5]

读取临床信息

#读取临床信息clin <- data.table::fread('Survival_SupplementalTable_S1_20171025_xena_sp',data.table = F)

筛选有用的临床资料，并将复杂的名字改成简单的名字

clin1 <-clin %>% select(sample,'cancer type abbreviation',gender,ajcc_pathologic_tumor_stage,OS,OS.time) %>% #选取这些列 rename(Type='cancer type abbreviation',Stage=ajcc_pathologic_tumor_stage,status=OS,time=OS.time)

读取基因注释文件（详见 TCGA数据下载与ID转换）

gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf.gz')#转换为数据框gtf <- as.data.frame(gtf)#保存save(gtf,file = 'gtf.Rdata')

去掉基因名小数点及后面的数字方便下一步转换

colnames(ALLdata)[1]<-'gene_id'ALLdata1<-separate(ALLdata,gene_id,into = c('gene_id'),sep='\\.') ALLdata1[1:5,1:5]

提取编码蛋白，与表达谱进行合并以转换基因名。（也可以取miRNA或者lncRNA等）

mRNA<-dplyr::filter(gtf,type=='gene',gene_biotype=='protein_coding')%>%#选择编码蛋白  select(gene_name,gene_id,gene_biotype)%>%#选择有用的三列  inner_join(ALLdata1,by ='gene_id')%>%#与表达谱合并  select(-gene_id,-gene_biotype)%>%  distinct(gene_name,.keep_all = T)mRNA[1:5,1:5]

将gene_name这一列作为行名

row.names(mRNA)<-mRNA[,1]mRNA<-mRNA[,-1]mRNA[1:5,1:5]

这个时候继续下去可能电脑hold不住了，先保存下文件

save(mRNA,file = 'pancancer_mRNA.Rdata')save(clin1,file = 'pancancer_clin.Rdata')

重新开始

rm(list = ls())load('pancancer_mRNA.Rdata')

将表达谱倒置，方便后续与临床资料的合并

mRNA<-as.data.frame(t(mRNA))mRNA[1:5,1:5]

接下来是设计Normal和Cancer的分组，根据TCGA数据ID的特性可区分Normal和Cancer TCGA差异分析及ggplot作图验证

group_list=ifelse(as.numeric(substr(rownames(mRNA),14,15)) < 10,'tumor','normal')design <- model.matrix(~0+factor(group_list))colnames(design)=levels(factor(group_list))rownames(design)=rownames(mRNA)head(design)

这里我用了之前做差异分析时分组的代码，应该可以写成更简洁的形式

给design数据数据变成数据框，添加一列ID

class(design)design<-as.data.frame(design)design$ID<-row.names(design)design[1:3,1:3]

将design分组中tumor这一列的1换成Tumor，0换成Normal

design$tumor[design$tumor==1] <- 'Tumor'#1换成Tumordesign$tumor[design$tumor==0] <- 'Normal'#0换成Normaldesign<-design[,-1]#去掉normal这一列colnames(design)[1]<-'Type'#cancer这一列列名改为Typehead(design)

合并数据

mRNA$ID<-row.names(mRNA)mRNA=inner_join(mRNA,design,by ='ID',copy=T)mRNA[1:5,1:5]

调整列的顺序便于观察

mRNA<-select(mRNA,ID,Type,everything())mRNA[1:5,1:5]

加载临床资料

load('pancancer_clin.Rdata')clin1[1:5,1:5]

为了合并，将列名sample换成ID，由于Type列名重复了，将Type换成Cancer

clin1<-rename(clin1,ID=sample,Cancer=Type)clin1[1:5,1:5]

合并

drawdata<-dplyr::inner_join(mRNA,clin1,by ='ID',copy=T)drawdata[1:5,1:5]

调整列名便于观察

drawdata<-select(drawdata,ID,Cancer,gender,Stage,status,time,everything())drawdata[1:5,1:10]

保存数据

save(drawdata,file = 'pancancer_drawdata.Rdata')

重新加载数据（以后就可以直接加载这个数据画图了，但是别忘了加载包）

rm(list = ls())load('pancancer_drawdata.Rdata')drawdata[1:5,1:10]

看看每个肿瘤和正常组织的个数

table(filter(drawdata,Type=='Tumor')$Cancer)table(filter(drawdata,Type=='Normal')$Cancer)

接下来就可以任意画图啦！

看一看单基因在泛癌中的表达

p<-ggplot(filter(drawdata,Type=='Tumor'),aes(x = Cancer, y =SPP1,color=Cancer))+ geom_boxplot()+ geom_jitter()+ theme_bw()+ theme(legend.position='none')
p

加上正常组织

library(ggpubr)p <- ggboxplot(drawdata, x = 'Cancer', y = 'SPP1',               fill = 'Type',legend=F,palette =c('#00AFBB', '#E7B800'))+  theme_bw()p

单个肿瘤差异表达

p <- ggboxplot(filter(drawdata,Cancer=='LIHC'), x = 'Type', y = 'SPP1', fill = 'Type',legend=F,palette =c('#00AFBB', '#E7B800'),bxp.errorbar=T)+ theme(legend.position='none')+ ylab(label = 'SPP1 expression')p

加上P值

my_comparisons <- list( c('Tumor', 'Normal') )p + stat_compare_means(comparisons = my_comparisons)

肿瘤分期表达

先看一下分期的情况

table(drawdata$Stage)

只取Stage I  II  III  IV期的数据，先整理分期情况

drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IVC' = 'Stage IV'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IVB' = 'Stage IV'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IVA' = 'Stage IV'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IIIC' = 'Stage III'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IIIB' = 'Stage III'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IIIA' = 'Stage III'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IIC' = 'Stage II'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IIB' = 'Stage II'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IIA' = 'Stage II'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IB' = 'Stage I'))drawdata$Stage<-str_replace_all(drawdata$Stage, c('Stage IA' = 'Stage I'))table(drawdata$Stage)

过滤数据只留下Stage I  II  III  IV期的数据

Stagedata<-filter(drawdata,Stage == 'Stage I' |Stage == 'Stage II' | Stage == 'Stage III' | Stage == 'Stage IV')table(Stagedata$Stage)

现在就可以画任意肿瘤任意基因分期的表达了

p <- ggboxplot(filter(Stagedata,Cancer=='LIHC'), x = 'Stage', y = 'SPP1',               fill = 'Stage',legend=F,bxp.errorbar=T)+  theme(legend.position='none')+  ylab(label = 'SPP1 expression')p

调整横坐标顺序

Stagedata$Stage<-factor(Stagedata$Stage,levels=c('Stage I','Stage II','Stage III','Stage IV'))p <- ggboxplot(filter(Stagedata,Cancer=='LIHC'), x = 'Stage', y = 'SPP1', fill = 'Stage',legend=F,bxp.errorbar=T)+ theme(legend.position='none')+ ylab(label = 'SPP1 expression')p

加上P值

my_comparisons <- list( c('Stage I', 'Stage II'),c('Stage I', 'Stage III'),c('Stage I', 'Stage IV'),c('Stage II', 'Stage III'),c('Stage II', 'Stage IV'),c('Stage III', 'Stage IV') )p + stat_compare_means(comparisons = my_comparisons,method = 't.test')

还可以看不同性别，该基因的差异表达，这里就不再举例子了，一般也用不到。

有了这个数据，我们还可以做任意基因在任意肿瘤中的生存分析

先加载包

library(survminer)library(survival)

看看Status和time这两列的情况

table(drawdata$status)table(drawdata$time)

把随访时间为0的，和正常的样本去掉

OSdata<-filter(drawdata,drawdata$time!=0 & drawdata$Type == 'Tumor' )table(OSdata$time)OSdata[1:5,1:10]

计算肝癌中SPP1表达的中位数

LIHC<-filter(OSdata,Cancer=='LIHC')med.exp<-median(LIHC$'SPP1')

根据中位数分成高表达和低表达两组

more.med.exp.index<-which(LIHC$'SPP1'>=med.exp)less.med.exp.index<-which(LIHC$'SPP1'<med.exp)LIHC$'SPP1'[more.med.exp.index]<-paste0('High expression (', length(more.med.exp.index),')')LIHC$'SPP1'[less.med.exp.index]<-paste0('Low expression (', length(less.med.exp.index),')')LIHC.fit<-survfit(Surv(time, status) ~ SPP1, data = LIHC)

画图

LIHC.survival<-ggsurvplot(LIHC.fit,                           data=LIHC,                           risk.table = T,                           risk.table.col = 'strata',                           risk.table.y.text = F,                           risk.table.title = 'Number at risk',                           pval = TRUE,                            #conf.int = TRUE,                           xlab = 'Time (days)',                           ggtheme = theme_light(),                           surv.median.line = 'hv',                           title='LIHC SPP1 Survival')
LIHC.survival

当然，我们还可以做任意一个分期，或者不同性别中任意基因的生存分析，操作方法都是一样的。

其实这个数据还可以做很多分析，包括任意两个基因在任意肿瘤中的相关性分析，以及在泛癌中的相关性分析等等。