DAY2-用TCGAbiolinks下载肝癌数据并做预处理1
DAY2-用TCGAbiolinks做生存分析2
正值春天，本该是踏青游玩的时刻，半个月来却阴雨连绵。今天终于见到了太阳，心情都变得愉悦了。今天继续接着上次做失败的生存分析，并打算学一学批量做生存分析。
# 先读取基因表达矩阵
setwd("E:\\MATH")
dataFilt_LIHC_final = read.csv("TCGA_LIHC_final.csv",row.names = 1,check.names=FALSE)
# 提取特定基因在癌症样本中的表达，选出的样本都是肿瘤样本
library(TCGAbiolinks)
samplesTP <- TCGAquery_SampleTypes(colnames(dataFilt_LIHC_final), typesample = c("TP"))
ABCB1 <- dataFilt_LIHC_final[c("ABCB1"),samplesTP]
## 修改样本名称（原是TCGA-DD-AAD5-01A-11R-A41C-07），但是clin.LIHC数据中的样本名称是12位，因为后期要把 clin.LIHC 数据和表达数据结合起来，需要将名字简化成TCGA-DD-AAD5
names(ABCB1) <- sapply(strsplit(names(ABCB1),'-'),function(x) paste0(x[1:3],collapse="-"))
#转置
ABCB1 <- t(ABCB1)
合并ABCB1基因与临床数据
# 合并ABCB1基因与临床数据
clin.LIHC <- GDCquery_clinic("TCGA-LIHC", "clinical")#重新下载临床数据
save(clin.LIHC,file = "clin.LIHC.Rdata")
load('clin.LIHC.Rdata')
clin.LIHC$"ABCB1" <- ABCB1[match(clin.LIHC$submitter_id,rownames(ABCB1)),]
df<-subset(clin.LIHC,select =c(submitter_id,vital_status,days_to_death,days_to_last_follow_up,ABCB1))
#去掉 NA
df <- df[!is.na(df$ABCB1),]
#根据ABCB1基因的表达情况进行分组，取基因的平均值，大于平均值的为H，小于为L
df$exp <- ''
df[df$ABCB1 >= mean(df$ABCB1),]$exp <- "H"
df[df$ABCB1 < mean(df$ABCB1),]$exp <- "L"
下面是df和最终版df
接下来分别用 TCGAanalyze_survival 函数和survival、survminer包进行生存分析，两者做比较。
首先使用TCGAanalyze_survival 函数进行单个基因表达情况对患者生存的影响。这个函数需要三个生存数据：患者的生存状态、死亡时间、随访时间。当存在失访的情况，就按照最后一次随访的时间来计算。
#用 TCGAanalyze_survival 函数进行生存分析
TCGAanalyze_survival(df,
clusterCol="exp",
risk.table = FALSE,
conf.int = FALSE,
color = c("Dark2"))
用 survival、survminer这两个包来进行生存分析.用status表示患者结局，1表示删失，2表示死亡。和使用TCGAanalyze_survival 函数不同的是，使用 survival 包来进行生存分析时，我们需要将 days_to_death 和 days_ to_last_follow_up 合并为一列。一般来说days_to_death 是空缺的， days_ to_last_follow_up 就会有数据，反之亦然，但是对比发现也有两列同时都有数据的，所以合并在一起是不是稍微有一点不妥呢？
df2 <- df
# 将status表示患者结局，1表示删失，2表示死亡
df2[df2$vital_status=='Dead',]$vital_status <- 2
df2[df2$vital_status=='Alive',]$vital_status <- 1
df2$vital_status <- as.numeric(df2$vital_status)
#将 days_to_death 和 days_ to_last_follow_up 合并为一列
df2$time <- df2$days_to_death
df2$time[which(is.na(df2$time))] <- df2$days_to_last_follow_up[which(is.na(df2$time))]
View(df2)
#进行生存分析
library(survival)
library(survminer)
# 建模
fit <- survfit(Surv(time, vital_status)~exp, data=df2) # 根据表达建模
# 显示P value
surv_pvalue(fit)$pval.txt
# 画图
ggsurvplot(fit,pval=TRUE)
#改版
ggsurvplot(
fit, # survfit object with calculated statistics.
data = df2, # data used to fit survival curves.
risk.table = TRUE, # show risk table.
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for point estimates of survival curves.
# survival estimates.
xlab = "Time", # customize X axis label.
ggtheme = theme_light(), # customize plot and risk table with a theme.
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.height = 0.25, # the height of the risk table
risk.table.y.text = FALSE,# show bars instead of names in text annotations in legend of risk table.
ncensor.plot = TRUE, # plot the number of censored subjects at time t
ncensor.plot.height = 0.25,
conf.int.style = "ribbon", # customize style of confidence intervals
surv.median.line = "hv", # add the median survival pointer.
)
改版后的图是不是更好看一些呢，最上面是生存曲线，配色都配好的，p值也标注好了，半数生存的时候对应的时间也标好了，95%置信区间也有了。中间一个表格是对应不同生存时间的患者数量，很直观的比较不同组别之间差异。下面一张图是不同分组的截断数据的数量，柱状图，简单明了。
有时候会在某一句代码那里卡住，需要琢磨很多天，或许在某一瞬间突然又可以做下去了，真是神奇。
参考：
