#' @name compSurv
#' @title Comparison of prognosis by Kaplan-Meier survival curve
#' @description This function calculates Kaplan-meier estimator and generate survival curves with log-rank test to detect prognostic difference among identified subtypes. If more than 2 subtypes are identified, pair-wise comparisons will be performed with an additional table printed on the survival curve.
#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.
#' @param surv.info A data.frame with rownames of observations and with at least two columns of `futime` for survival time and `fustat` for survival status (0: censoring; 1: event)
#' @param clust.col A string vector storing colors for each subtype.
#' @param p.adjust.method A string value for indicating method for adjusting p values (see \link[stats]{p.adjust}). Allowed values include one of c(`holm`, `hochberg`, `hommel`, `bonferroni`, `BH`, `BY`, `fdr`, `none`); "BH" by default.
#' @param fig.name A string value to indicate the output path for storing the kaplan-meier curve.
#' @param fig.path A string value to indicate the name of the kaplan-meier curve.
#' @param convt.time A string value to indicate how to convert the survival time; value of `d` for days, `m` for months and `y` for years; "d" by default.
#' @param surv.median.line A string value for drawing a horizontal/vertical line at median survival. Allowed values include one of c(`none`, `hv`, `h`, `v`). v: vertical, h:horizontal; "none" by default.
#' @param xyrs.est An integer vector to estimate probability of surviving beyond a certain number (x) of years (Estimating x-year survival); NULL by default.
#' @param surv.cut A numeric value to indicate the x-axis cutoff for showing the maximal survival time. NULL by default (show 0-maximum survival time range).
#' @return A figure of multi-omics Kaplan-Meier curve (.pdf) and a list with the following components:
#'
#' \code{fitd} an object returned by \link[survival]{survdiff}.
#'
#' \code{fid} an object returned by \link[survival]{survfit}.
#'
#' \code{xyrs.est} x-year probability of survival and the associated lower and upper bounds of the 95% confidence interval are also displayed if argument of `xyrs.est` was set by users.
#'
#' \code{overall.p} a nominal p.value calculated by Kaplain-Meier estimator with log-rank test.
#'
#' \code{pairwise.p} an object of class "pairwise.htest" which is a list containing the p values (see \link[survminer]{pairwise_survdiff}); (\emph{only returned when more than 2 subtypes are identified}).
#' @import survival
#' @import survminer
#' @import ggplot2
#' @importFrom grDevices pdf dev.off pdf.options
#' @importFrom ggpp geom_table
#' @importFrom tibble tibble
#' @export
#' @examples # There is no example and please refer to vignette.
compSurv <- function(moic.res = NULL,
surv.info = NULL,
convt.time = "d",
surv.cut = NULL,
xyrs.est = NULL,
clust.col = c("#2EC4B6","#E71D36","#FF9F1C","#BDD5EA","#FFA5AB","#011627","#023E8A","#9D4EDD"),
p.adjust.method = "BH",
surv.median.line = "none",
fig.name = NULL,
fig.path = getwd()){
if(!all(is.element(c("futime","fustat"), colnames(surv.info)))) {
stop("fail to find variables of futime and fustat.")
}
if(!all(is.element(convt.time, c("d","m","y")))) {
stop("unsupported time conversion. Allowed values contain c('d', 'm', 'y').")
}
# get common samples
comsam <- intersect(rownames(surv.info),rownames(moic.res$clust.res))
mosurv.res <- cbind.data.frame(surv.info[comsam,c("futime","fustat")],
moic.res$clust.res[comsam, "clust", drop = FALSE])
message(paste0("--a total of ",length(comsam), " samples are identified."))
# remove missing data if possible
if(sum(c(is.na(mosurv.res$futime), is.na(mosurv.res$fustat))) > 0) {
message("--removed missing values.")
mosurv.res <- as.data.frame(na.omit(mosurv.res))
message(paste0("--leaving ",nrow(mosurv.res), " observations."))
}
# check time unit
if(max(mosurv.res$futime) < 365) {
warning("it seems the 'futime' might not in [day] unit, please make sure you provide the correct survival information.")
}
# create new variable of Subtype
mosurv.res$Subtype <- paste0("CS", mosurv.res$clust)
mosurv.res <- mosurv.res[order(mosurv.res$Subtype),]
# estimate x-year survival probability
if(!is.null(xyrs.est)) {
if(max(xyrs.est) * 365 < max(mosurv.res$futime)) {
xyrs <- summary(survfit(Surv(futime, fustat) ~ Subtype, data = mosurv.res), times = xyrs.est * 365)
} else {
stop("the maximal survival time is less than the time point you want to estimate!")
}
} else {
xyrs <- "[Not Available]: argument of xyrs.est was not specified."
}
# convert survival time
mosurv.res$futime <- switch(convt.time,
"d" = mosurv.res$futime,
"m" = round(mosurv.res$futime/30.5,4),
"y" = round(mosurv.res$futime/365,4))
date.lab <- switch(convt.time,
"d" = "Days",
"m" = "Months",
"y" = "Years")
# truncate survival time
if(date.lab == "Days" & max(mosurv.res$futime) > 3650) {
#xlim = c(0, 3650)
brk = 365
}
if(date.lab == "Days" & max(mosurv.res$futime) <= 3650) {
#xlim = c(0, max(mosurv.res$futime))
brk = floor(max(mosurv.res$futime)/10)
}
if(date.lab == "Months" & max(mosurv.res$futime) > 120) {
#xlim = c(0, 120)
brk = 12
}
if(date.lab == "Months" & max(mosurv.res$futime) <= 120) {
#xlim = c(0, max(mosurv.res$futime))
brk = floor(max(mosurv.res$futime)/10)
}
if(date.lab == "Years" & max(mosurv.res$futime) > 10) {
#xlim = c(0, 10)
brk = 1
}
if(date.lab == "Years" & max(mosurv.res$futime) <= 10) {
#xlim = c(0, max(mosurv.res$futime))
brk = 1
}
if(is.null(surv.cut)) {
xlim = c(0, max(mosurv.res$futime))
} else {
message(paste0("--cut survival curve up to ",surv.cut," ",tolower(date.lab)))
xlim = c(0, surv.cut)
}
# if(!is.null(surv.cut)) {
# xlim = c(0, surv.cut)
# brk = surv.cut/10
# } else {message("--cut survival curve up to 10 years.")}
n.moic <- length(unique(mosurv.res$Subtype))
# basic survival analysis
fitd <- survdiff(Surv(futime, fustat) ~ Subtype,
data = mosurv.res,
na.action = na.exclude)
p.val <- 1 - pchisq(fitd$chisq, length(fitd$n) - 1)
fit <- survfit(Surv(futime, fustat)~ Subtype,
data = mosurv.res,
type = "kaplan-meier",
error = "greenwood",
conf.type = "plain",
na.action = na.exclude)
# hack strata for better survival curve
names(fit$strata) <- gsub("Subtype=", "", names(fit$strata))
# kaplan-meier curve
p <- suppressWarnings(ggsurvplot(fit = fit,
conf.int = FALSE,
risk.table = TRUE,
risk.table.col = "strata",
palette = clust.col[1:n.moic],
data = mosurv.res,
size = 1,
xlim = xlim,
break.time.by = brk,
legend.title = "",
surv.median.line = surv.median.line,
xlab = paste0("Time (",date.lab,")"),
ylab = "Survival probability (%)",
risk.table.y.text = FALSE))
# make survival time as percentage
p$plot <- quiet(p$plot + scale_y_continuous(breaks = seq(0, 1, 0.25), labels = seq(0,100,25)))
if(n.moic > 2) {
# add nominal pvalue for log-rank test
p.lab <- paste0("Overall P",
ifelse(p.val < 0.001, " < 0.001",
paste0(" = ",round(p.val, 3))))
p$plot <- p$plot + annotate("text",
x = 0, y = 0.55,
hjust = 0,
fontface = 4,
label = p.lab)
# calculate pair-wise survival comparison
ps <- pairwise_survdiff(Surv(futime, fustat)~ Subtype,
data = mosurv.res,
p.adjust.method = p.adjust.method)
# add pair-wise comparison table
# options(stringsAsFactors = FALSE)
addTab <- as.data.frame(as.matrix(ifelse(round(ps$p.value, 3) < 0.001, "<0.001",
round(ps$p.value, 3))))
addTab[is.na(addTab)] <- "-"
# options(stringsAsFactors = TRUE)
df <- tibble(x = 0, y = 0, tb = list(addTab))
p$plot <- p$plot + geom_table(data = df, aes(x = x, y = y, label = tb), table.rownames = TRUE)
} else {
# add nominal pvalue for log-rank test
p.lab <- paste0("P",
ifelse(p.val < 0.001, " < 0.001",
paste0(" = ",round(p.val, 3))))
p$plot <- p$plot + annotate("text",
x = 0, y = 0.55,
hjust = 0,
fontface = 4,
label = p.lab)
}
# output curve to pdf
if(!is.null(fig.name)) {
outFile <- file.path(fig.path,paste0(fig.name,".pdf"))
} else {
outFile <- file.path(fig.path,paste0("km_curve_",moic.res$mo.method,".pdf"))
}
pdf.options(reset = TRUE, onefile = FALSE)
pdf(outFile, width = 6, height = 7)
# ggsave(outFile, width = 6, height = 7)
print(p)
dev.off()
# output curve to screen
print(p)
if(n.moic > 2) {
return(list(fitd = fitd, fit = fit, xyrs.est = xyrs, overall.p = p.val, pairwise.p = ps))
} else {
return(list(fitd = fitd, fit = fit, xyrs.est = xyrs, overall.p = p.val))
}
}
Datasets
Models
