TCGA-mRNA的生存分析

引言：本文主要记录mRNA生存分析过程中遇到的问题以及具体过程，生存分析用的是最常见的KM法Kaplan-Meier。

工具：
TCGA的数据通过library(TCGAbiolinks)下载：癌症的clinical数据、miRNA/mRNA的表达量数据
生存分析由library(sirvival)完成

分析过程中问题：

参考脚本：http://www.zxzyl.com/archives/1063

1.结果图中不显示加号，不显示删失数据

解决：
clincial文件注释：https://www.cnblogs.com/emanlee/p/7635951.html
其中有三列：
day_to_death:距离死亡的时间
day_to_last_follow_up：最后一个随访时间到开始时间(Time interval from the date of last followup to the date of initial pathologic diagnosis, represented as a calculated number of days.)
vital_status是生存状态，alive和dead
alive表示删失数据，dead表示发生了终点事件。

days_to_death和days_to_last_follow_up.png

图中为day_to_death(左边)和day_to_last_follow_up，可以看出两列基本互补，对于dead样本基本对应的有days_to_death,而基本都没有days_to_last_follow_up,对于alive的数据则是没有days_to_death数据，有days_to_last_follow_up。

也有的人是将NA变成0，然后两列相加作为最终的OS.这里发现有两列都为NA的情况，也有两列都有数值的情况。

所以，最后采用的是：
对于dead样本，overall survival采用day_to_death
对于alive样本，overall survival采用day_to_last_follow_up

df$os<-ifelse(df$vital_status=='Alive',df$days_to_last_follow_up,df$days_to_death)

2.高表达和低表达mRNA的对生存的影响
整合所有癌症样本中特定gene的表达量，大于表达量均值的为高表达，低于表达量均值的为低表达，根据样本号和clinical文件中的submitter_id比对，找出这些样本的生存时间OS

R-Script-mRNA

library(SummarizedExperiment)
library(TCGAbiolinks)
library(survival)
library(survminer)
TCGAbiolinks:::getProjectSummary("TCGA-LIHC") 
clinical <- GDCquery_clinic(project = "TCGA-LIHC", type = "clinical")
query <- GDCquery(project = "TCGA-LIHC", 
                  experimental.strategy = "RNA-Seq",
                  data.category = "Transcriptome Profiling", 
                  data.type = "Gene Expression Quantification",
                  workflow.type = "HTSeq - FPKM")
GDCdownload(query)
Rnaseq <- GDCprepare(query)
fpkm_matrix <- assay(Rnaseq) #【fig：gene矩阵】
#---------第一部分结束，下载了mRNA测序数据以及clinical数据---------------------
samplesTP <- TCGAquery_SampleTypes(barcode = colnames(fpkm_matrix),typesample = c("TP"))
#这里可以去查TCGAquery_SampleTypes()的参数说明，这一步找出所有的癌症样本
gene_exp <- fpkm_matrix[c("ENSGxxxxx"),samplesTP]
#传入一个gene的ENSEMBL_ID,从所有gene的多个样本中筛选出这个基因的特定样本，给定选取特定行和列的参数。
names(gene_exp) <-  sapply(strsplit(names(gene_exp),'-'),function(x) paste(x[1:3],collapse="-"))
#把RNA-seq中的列名保留前三个，因为clinical文件中的submitter id只有前三段
clinical$GENE <- gene_exp[clinical$submitter_id]
#给clinical文件新加一列GENE，就是我们研究的基因表达量
#-----------------第二部分结束------整合完最终需要的文件-------------
df<-subset(clinical,select=c(submitter_id,vital_status,days_to_death,days_to_last_follow_up,GENE))
df$os<-ifelse(df$vital_status=='Alive',df$days_to_last_follow_up,df$days_to_death)
#alive的样本采用days_to_last_follow_up
df <- df[!is.na(df$GENE),]#去掉表达量为0的样本
df$exp=''
df[df$GENE >= mean(df$GENE),]$exp <- "High"
df[df$GENE <  mean(df$GENE),]$exp <- "Low"
df[df$vital_status=='Dead',]$vital_status <- 2
df[df$vital_status=='Alive',]$vital_status <- 1
df$vital_status <- as.numeric(df$vital_status)
#--------------第三部分结束----对mRNA的高低表达量进行确定------
fit <- survfit(Surv(os, vital_status)~exp, data=df) # 根据表达建模
head(summary(fit))
ggsurvplot(fit,pval=TRUE)

gene矩阵.png

结果图(就是随手画了个这么没差异的gene)：

Rplot01.png

a<-ggsurvplot(fit,legend.title = "Expression",palette = c('red','black'), 
           pval = TRUE,
           risk.table = TRUE,
           tables.height = 0.2,
           tables.theme = theme_cleantable(),
           xlab='Time(days)'
)
ggsave('/home/zhang/xxxx/',print(a)) #输出图片

re.png