DAY3-跟着另一篇文章做生存分析
总体思路：TCGAbiolinks包下载肝癌mRNA的FPKM类型数据及临床数据，用survival包做生存分析。
1.在clinical文件中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。大多数把这两列相加作为最终的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。
#数据下载并合成矩阵
setwd("E:\\MATH")
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矩阵】
save(fpkm_matrix,file = "fpkm_matrix.Rdata")#保存文件
以下是下载后数据的格式
clinical.png
Rnaseq.png
fpkm_matrix.png
#整合clinical和基因表达矩阵文件
#这里可以去查TCGAquery_SampleTypes()的参数说明，这一步找出所有的癌症样本
samplesTP <- TCGAquery_SampleTypes(barcode = colnames(fpkm_matrix),typesample = c("TP"))
#传入一个gene的ENSEMBL_ID,从所有gene的多个样本中筛选出这个基因的特定样本，给定选取特定行和列的参数。我随便选了一个
gene_exp <- fpkm_matrix[c("ENSG00000000003"),samplesTP]
#把RNA-seq中的列名保留前三个，因为clinical文件中的submitter id只有前三段
names(gene_exp) <- sapply(strsplit(names(gene_exp),'-'),function(x) paste(x[1:3],collapse="-"))
#给clinical文件新加一列GENE，就是我们研究的基因表达量
clinical$GENE <- gene_exp[clinical$submitter_id]
#在clinical中提取我们想要的列
df<-subset(clinical,select=c(submitter_id,vital_status,days_to_death,days_to_last_follow_up,GENE))
#alive的样本采用days_to_last_follow_up
df$os<-ifelse(df$vital_status=='Alive',df$days_to_last_follow_up,df$days_to_death)
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)
最终的df.png
#画图
fit <- survfit(Surv(os, vital_status)~exp, data=df) # 根据表达建模
head(summary(fit))
ggsurvplot(fit,pval=TRUE)
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)'
)
生存曲线.png
（放很多图是怕忘记这种结构，跟别的对比的时候想不起来）
觉得很多时间都花费在下载数据，对于生存分析也是似懂非懂，多实践就可以总结出自己的思想了。又好像有点明白昨天为什么做不出生存曲线了，但是又抓不住那一闪而过的想法，可能是夜深了，大脑不运转了。而且今天的图感觉可以做得更美。还是有很多问题，虽然是下载数据合并成矩阵了，是否还应该对数据进行下一步处理，要是我想知道这个基因编号，难道我要到ensembl去查才知道吗？这是下载FPKM数据，要是下载HTSeq - Counts数据呢又会怎样？
