source(file = "/home/longzhilin/Analysis_Code/SurvivalAnalysis/Cox.function.R")
RCC.icb.analysis <- function(signature.list, expresionMatrix, clincal.info){
library(ggpubr)
signature.list.score <- lapply(names(signature.list), function(y){
genes <- signature.list[[y]]
programName <- y
index <- match(genes, rownames(expresionMatrix))
if(length(na.omit(index)) < length(index)){
cat(y, " lack ", genes[which(is.na(index))], " gene\n")
}
signature.score <- colMeans(expresionMatrix[na.omit(index),])
ICB.treatment.idx <- which(clincal.info$Arm == "NIVOLUMAB")
signature.score.icb <- signature.score[ICB.treatment.idx]
clincal.info.icb <- clincal.info[ICB.treatment.idx,]
plot.data.icb <- data.frame(score = signature.score.icb, Cohort = clincal.info.icb$Cohort, Benefit = clincal.info.icb$Benefit)
p <- ggboxplot(plot.data.icb, x = "Cohort", y = "score", color = "Benefit", xlab = "", ylab = "Signature score", title = paste0("ICB threapy", "-", programName), add = "jitter") + stat_compare_means(aes(group = Benefit), label = "p.format")
print(p)
index <- which(plot.data.icb$Benefit=="ICB")
plot.data.icb.CB.VS.ICB <- plot.data.icb[-index,]
p <- ggboxplot(plot.data.icb.CB.VS.ICB, x = "Cohort", y = "score", color = "Benefit", xlab = "", ylab = "Signature score", title = paste0("ICB threapy", "-", programName), add = "jitter") + stat_compare_means(aes(group = Benefit), label = "p.format")
print(p)
mTOR.treatment.idx <- which(clincal.info$Arm == "EVEROLIMUS")
signature.score.mTOR <- signature.score[mTOR.treatment.idx]
clincal.info.mTOR <- clincal.info[mTOR.treatment.idx,]
plot.data.mTOR <- data.frame(score = signature.score.mTOR, Cohort = clincal.info.mTOR$Cohort, Benefit = clincal.info.mTOR$Benefit)
p <- ggboxplot(plot.data.mTOR, x = "Cohort", y = "score", color = "Benefit", xlab = "", ylab = "Signature score", title = paste0("mTOR threapy", "-", programName), add = "jitter") + stat_compare_means(aes(group = Benefit), label = "p.format")
print(p)
index <- which(plot.data.mTOR$Benefit=="ICB")
plot.data.mTOR.CB.VS.ICB <- plot.data.mTOR[-index,]
p <- ggboxplot(plot.data.mTOR.CB.VS.ICB, x = "Cohort", y = "score", color = "Benefit", xlab = "", ylab = "Signature score", title = paste0("mTOR threapy", "-", programName), add = "jitter") + stat_compare_means(aes(group = Benefit), label = "p.format")
print(p)
#survival plot
#Based on median grouping
survival.res <- sapply(c("CM-009", "CM-010", "CM-025"), function(x){
index <- which(clincal.info.icb$Cohort == x)
clincal.info.icb.cohort <- clincal.info.icb[index,]
signature.score.icb.cohort <- signature.score.icb[index]
mid.score.icb <- median(signature.score.icb.cohort)
icb.label <- rep("High score", length(signature.score.icb.cohort))
icb.label[which(signature.score.icb.cohort<=mid.score.icb)] <- "Low score"
if(length(table(icb.label)) <= 1){
cox.res <- "NA"
}else{
if(x == "CM-025"){
# different treatment
icb.clincal.OS <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, event = clincal.info.icb.cohort$OS_CNSR, time = clincal.info.icb.cohort$OS, sample.label = icb.label)
p1 <- plot.surv(icb.clincal.OS, risk.table = T, HR = T, ylab = "Overall Survival", main = paste0("ICB ", programName, " ", x), surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
icb.clincal.PFS <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, event = clincal.info.icb.cohort$PFS_CNSR, time = clincal.info.icb.cohort$PFS, sample.label = icb.label)
p1 <- plot.surv(icb.clincal.PFS, risk.table = T, HR = T, ylab = "Progression-Free Survival", main = paste0("ICB ", programName, " ", x), surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
cox.clincal.data <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, time = clincal.info.icb.cohort$OS, event = clincal.info.icb.cohort$OS_CNSR, sample.label = icb.label, Sex = clincal.info.icb.cohort$Sex, Age = clincal.info.icb.cohort$Age, Prior_Therapies = clincal.info.icb.cohort$Number_of_Prior_Therapies, Metastasis = clincal.info.icb.cohort$Tumor_Sample_Primary_or_Metastasis)
cox.clincal.data <- na.omit(cox.clincal.data)
icb.cox.OS <- Cox.function(time = cox.clincal.data$time, event = cox.clincal.data$event, clinical.data = cox.clincal.data)
cox.clincal.data <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, time = clincal.info.icb.cohort$PFS, event = clincal.info.icb.cohort$PFS_CNSR, sample.label = icb.label, Sex = clincal.info.icb.cohort$Sex, Age = clincal.info.icb.cohort$Age, Prior_Therapies = clincal.info.icb.cohort$Number_of_Prior_Therapies, Metastasis = clincal.info.icb.cohort$Tumor_Sample_Primary_or_Metastasis)
cox.clincal.data <- na.omit(cox.clincal.data)
icb.cox.PFS <- Cox.function(time = cox.clincal.data$time, event = cox.clincal.data$event, clinical.data = cox.clincal.data)
####mTOR
mid.score.mTOR <- median(signature.score.mTOR)
mTOR.label <- rep("High score", length(signature.score.mTOR))
mTOR.label[which(signature.score.mTOR<=mid.score.mTOR)] <- "Low score"
mTOR.clincal.OS <- data.frame(Patient_ID = clincal.info.mTOR$RNA_ID, event = clincal.info.mTOR$OS_CNSR, time = clincal.info.mTOR$OS, sample.label = mTOR.label)
p1 <- plot.surv(mTOR.clincal.OS, risk.table = T, HR = T, ylab = "Overall Survival", main = paste0("mTOR ", programName, " ", x), surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
mTOR.clincal.PFS <- data.frame(Patient_ID = clincal.info.mTOR$RNA_ID, event = clincal.info.mTOR$PFS_CNSR, time = clincal.info.mTOR$PFS, sample.label = mTOR.label)
p1 <- plot.surv(mTOR.clincal.PFS, risk.table = T, HR = T, ylab = "Progression-Free Survival", main = paste0("mTOR ", programName, " ", x), surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
cox.clincal.data <- data.frame(Patient_ID = clincal.info.mTOR$RNA_ID, time = clincal.info.mTOR$OS, event = clincal.info.mTOR$OS_CNSR, sample.label = mTOR.label, Sex = clincal.info.mTOR$Sex, Age = clincal.info.mTOR$Age, Prior_Therapies = clincal.info.mTOR$Number_of_Prior_Therapies, Metastasis = clincal.info.mTOR$Tumor_Sample_Primary_or_Metastasis)
cox.clincal.data <- na.omit(cox.clincal.data)
mTOR.cox.OS <- Cox.function(time = cox.clincal.data$time, event = cox.clincal.data$event, clinical.data = cox.clincal.data)
cox.clincal.data <- data.frame(Patient_ID = clincal.info.mTOR$RNA_ID, time = clincal.info.mTOR$PFS, event = clincal.info.mTOR$PFS_CNSR, sample.label = mTOR.label, Sex = clincal.info.mTOR$Sex, Age = clincal.info.mTOR$Age, Prior_Therapies = clincal.info.mTOR$Number_of_Prior_Therapies, Metastasis = clincal.info.mTOR$Tumor_Sample_Primary_or_Metastasis)
cox.clincal.data <- na.omit(cox.clincal.data)
mTOR.cox.PFS <- Cox.function(time = cox.clincal.data$time, event = cox.clincal.data$event, clinical.data = cox.clincal.data)
cox.res <- list(icb.OS = icb.cox.OS, icb.PFS = icb.cox.PFS, mTOR.OS = mTOR.cox.OS, mTOR.PFS = mTOR.cox.PFS)
}else{
icb.clincal.OS <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, event = clincal.info.icb.cohort$OS_CNSR, time = clincal.info.icb.cohort$OS, sample.label = icb.label)
p1 <- plot.surv(icb.clincal.OS, risk.table = T, HR = T, ylab = "Overall Survival", main = paste0(programName, " ", x), surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
icb.clincal.PFS <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, event = clincal.info.icb.cohort$PFS_CNSR, time = clincal.info.icb.cohort$PFS, sample.label = icb.label)
p1 <- plot.surv(icb.clincal.PFS, risk.table = T, HR = T, ylab = "Progression-Free Survival", main = paste0(programName, " ", x), surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
cox.clincal.data <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, time = clincal.info.icb.cohort$OS, event = clincal.info.icb.cohort$OS_CNSR, sample.label = icb.label, Sex = clincal.info.icb.cohort$Sex, Age = clincal.info.icb.cohort$Age)
cox.clincal.data <- na.omit(cox.clincal.data)
icb.cox.OS <- Cox.function(time = cox.clincal.data$time, event = cox.clincal.data$event, clinical.data = cox.clincal.data)
cox.clincal.data <- data.frame(Patient_ID = clincal.info.icb.cohort$RNA_ID, time = clincal.info.icb.cohort$PFS, event = clincal.info.icb.cohort$PFS_CNSR, sample.label = icb.label, Sex = clincal.info.icb.cohort$Sex, Age = clincal.info.icb.cohort$Age)
cox.clincal.data <- na.omit(cox.clincal.data)
icb.cox.PFS <- Cox.function(time = cox.clincal.data$time, event = cox.clincal.data$event, clinical.data = cox.clincal.data)
cox.res <- list(icb.OS = icb.cox.OS, icb.PFS = icb.cox.PFS)
}
}
return(cox.res)
})
})
names(signature.list.score) <- names(signature.list)
return(signature.list.score)
}
plot.surv <- function(clinical.data, upper.time = NULL, xscale = 1, xlab = "Time", median.time = TRUE,
surv.median.line = "none", HR = FALSE, risk.table = TRUE, pval = TRUE,
conf.int = FALSE, main = NULL, ylab = "Survival probability", colors = c("#D95319", "#3B6793","#EA4335","#4285F4","#34A853","#000000")) {
#Load related R packages
require(survival)
require(survminer)
require(RColorBrewer)
require(gridExtra)
#Determine the unit of event type and time
# survival.event <- survival.event[1];
# unit.xlabel <- unit.xlabel[1];
#If upper.time is set, the samples whose survival time exceeds upper.time will be removed
if (!is.null(upper.time)) clinical.data <- clinical.data[clinical.data$time <= upper.time,]
#set color
if (!is.factor(clinical.data$sample.label))
clinical.data$sample.label <- as.factor(clinical.data$sample.label)
t.name <- levels(clinical.data$sample.label)
if (length(t.name)> 6) stop("Sample grouping>6, exceeding the function acceptance range")
t.col <- colors[1:length(t.name)]
# Construct a living object
km.curves <- survfit(Surv(time, event)~sample.label, data=clinical.data)
# Calculate HR value and 95% CI
if (length(t.name) == 2) {
if (HR) {
cox.obj <- coxph(Surv(time, event)~sample.label, data=clinical.data)
tmp <- summary(cox.obj)
HRs <- round(tmp$coefficients[ ,2], digits = 2)
HR.confint.lower <- round(tmp$conf.int[,"lower .95"], 2)
HR.confint.upper <- round(tmp$conf.int[,"upper .95"], 2)
HRs <- paste0(HRs, " (", HR.confint.lower, "-", HR.confint.upper, ")")
}
}
# Construct the legend display text in the survival image
legend.content <- substr(names(km.curves$strata), start = 14, stop = 1000)
# x-axis scale unit conversion
if (is.numeric(xscale) | (xscale %in% c("d_m", "d_y", "m_d", "m_y", "y_d", "y_m"))) {
xscale = xscale
} else {
stop('xscale should be numeric or one of c("d_m", "d_y", "m_d", "m_y", "y_d", "y_m").')
}
# Implicit function: conversion of survival time unit
.format_xticklabels <- function(labels, xscale){
# 1 year = 365.25 days
# 1 month = 365.25/12 = 30.4375 days
if (is.numeric(xscale)) xtrans <- 1/xscale
else
xtrans <- switch(xscale,
d_m = 12/365.25,
d_y = 1/365.25,
m_d = 365.25/12,
m_y = 1/12,
y_d = 365.25,
y_m = 12,
1
)
round(labels*xtrans,2)
}
# Add the median survival time and its 95% CI in the figure and place it in the subtitle position
subtitle <- NULL
if (median.time) {
if (is.numeric(xscale)) {
median.km.obj = km.curves
} else if (xscale %in% c("d_m", "d_y", "m_d", "m_y", "y_d", "y_m")) {
clinical.data$time <- .format_xticklabels(labels = clinical.data$time, xscale = xscale)
median.km.obj <- survfit(Surv(time, event)~sample.label, data=clinical.data)
}
survival.time.info <- NULL
survival.time.info <- rbind(survival.time.info, summary(median.km.obj)$table)
median.survival <- round(survival.time.info[!duplicated(survival.time.info[,7:9]),7:9], digits = 2) # 注意：这里取得的置信区间上界可能为NA
if (length(levels(clinical.data$sample.label)) == 1) {
tmp1 <- levels(clinical.data$sample.label)
} else {
tmp1 <- do.call(rbind,strsplit(rownames(summary(median.km.obj)$table), split = "="))[,2]
}
tmp2 <- paste(median.survival[,1], "(", median.survival[,2], "-", median.survival[,3], ")")
subtitle <- paste(tmp1, tmp2, sep = ":", collapse = "\n")
}
# ggsurvplot
ggsurv <- ggsurvplot(km.curves, # survfit object with calculated statistics.
data = clinical.data, # data used to fit survival curves.
palette = t.col,
risk.table = risk.table, # show risk table.
pval = pval, # show p-value of log-rank test.
surv.median.line = surv.median.line, # add the median survival pointer.
title = main, #main title
subtitle = subtitle, #sub title
font.main = 15,
xlab = xlab, # customize X axis label.
ylab = ylab, # customize Y axis label
xscale = xscale,
#legend
legend.title = "",
legend.labs = legend.content,
legend = c(0.8,0.9),
font.legend = 9,
#risk table
tables.theme = theme_cleantable(),
risk.table.title = "No. at risk:",
risk.table.y.text.col = T,
risk.table.y.text = FALSE,
tables.height = 0.15,
risk.table.fontsize = 3
)
if (length(t.name) == 2) {
if (HR)
ggsurv$plot <- ggsurv$plot + ggplot2::annotate("text", x = max(km.curves$time)/12,
y = 0.15, size = 5, label = paste("HR=", HRs))
}
ggsurv$plot <- ggsurv$plot + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(size = 10),
plot.margin = unit(c(5.5, 5.5, 5.5, 50), "points"))
ggsurv$table <- ggsurv$table + theme(plot.title = element_text(hjust = -0.04),
plot.margin = unit(c(5.5, 5.5, 5.5, 50), "points"))
if(length(t.name) > 2) {
# pairwise: log rank P value
res <- pairwise_survdiff(Surv(time, event)~sample.label, data=clinical.data);
pairwise.pvalue <- round(res$p.value, digits = 4);
pairwise.pvalue[which(pairwise.pvalue < 0.0001)] <- "<0.0001";
pairwise.pvalue[is.na(pairwise.pvalue)] <- "-"
# add table
tt <- ttheme_minimal(core = list(fg_params = list(col = "black"),bg_params = list(fill = NA, col = "black")),
colhead = list(fg_params = list(col = NA),bg_params = list(fill = t.col, col = "black")),
rowhead = list(fg_params = list(col = NA, hjust = 1),bg_params = list(fill = c("white",t.col[-1]), col = "black"))
)
pairwise.table <- tableGrob(pairwise.pvalue, theme = tt)
ggsurv <- ggarrange(ggarrange(ggsurv$plot, ggsurv$table, nrow=2, heights=c(2,0.5)),
pairwise.table, nrow=2, heights = c(2,0.5),
labels = c("","p from pairwise comparisons"),
hjust = 0, font.label = list(size = 15, face = "plain"))
}
ggsurv
}
Datasets
Models
