引言：本文主要記錄mRNA生存分析過程中遇到的問題以及具體過程，生存分析用的是最常見的KM法Kaplan-Meier。
工具：
TCGA的數(shù)據(jù)通過library(TCGAbiolinks)下載：癌癥的clinical數(shù)據(jù)、miRNA/mRNA的表達(dá)量數(shù)據(jù)
生存分析由library(sirvival)完成
分析過程中問題：
參考腳本：http://www.zxzyl.com/archives/1063
1.結(jié)果圖中不顯示加號(hào)，不顯示刪失數(shù)據(jù)
解決：
clincial文件注釋：https://www.cnblogs.com/emanlee/p/7635951.html
其中有三列：
day_to_death:距離死亡的時(shí)間
day_to_last_follow_up：最后一個(gè)隨訪時(shí)間到開始時(shí)間(Time interval from the date of last followup to the date of initial pathologic diagnosis, represented as a calculated number of days.)
vital_status是生存狀態(tài)，alive和dead
alive表示刪失數(shù)據(jù)，dead表示發(fā)生了終點(diǎn)事件。
圖中為day_to_death(左邊)和day_to_last_follow_up，可以看出兩列基本互補(bǔ)，對(duì)于dead樣本基本對(duì)應(yīng)的有days_to_death,而基本都沒有days_to_last_follow_up,對(duì)于alive的數(shù)據(jù)則是沒有days_to_death數(shù)據(jù)，有days_to_last_follow_up。
也有的人是將NA變成0，然后兩列相加作為最終的OS.這里發(fā)現(xiàn)有兩列都為NA的情況，也有兩列都有數(shù)值的情況。
所以，最后采用的是：
對(duì)于dead樣本，overall survival采用day_to_death
對(duì)于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.高表達(dá)和低表達(dá)mRNA的對(duì)生存的影響
整合所有癌癥樣本中特定gene的表達(dá)量，大于表達(dá)量均值的為高表達(dá)，低于表達(dá)量均值的為低表達(dá)，根據(jù)樣本號(hào)和clinical文件中的submitter_id比對(duì)，找出這些樣本的生存時(shí)間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矩陣】
#---------第一部分結(jié)束，下載了mRNA測(cè)序數(shù)據(jù)以及clinical數(shù)據(jù)---------------------
samplesTP <- TCGAquery_SampleTypes(barcode = colnames(fpkm_matrix),typesample = c("TP"))
#這里可以去查TCGAquery_SampleTypes()的參數(shù)說(shuō)明，這一步找出所有的癌癥樣本
gene_exp <- fpkm_matrix[c("ENSGxxxxx"),samplesTP]
#傳入一個(gè)gene的ENSEMBL_ID,從所有g(shù)ene的多個(gè)樣本中篩選出這個(gè)基因的特定樣本，給定選取特定行和列的參數(shù)。
names(gene_exp) <- sapply(strsplit(names(gene_exp),'-'),function(x) paste(x[1:3],collapse="-"))
#把RNA-seq中的列名保留前三個(gè)，因?yàn)閏linical文件中的submitter id只有前三段
clinical$GENE <- gene_exp[clinical$submitter_id]
#給clinical文件新加一列GENE，就是我們研究的基因表達(dá)量
#-----------------第二部分結(jié)束------整合完最終需要的文件-------------
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),]#去掉表達(dá)量為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)
#--------------第三部分結(jié)束----對(duì)mRNA的高低表達(dá)量進(jìn)行確定------
fit <- survfit(Surv(os, vital_status)~exp, data=df) # 根據(jù)表達(dá)建模
head(summary(fit))
ggsurvplot(fit,pval=TRUE)
結(jié)果圖(就是隨手畫了個(gè)這么沒差異的gene)：
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)) #輸出圖片
