R/compSurv.R

In xlucpu/MOVICS: Multi-Omics integration and VIsualization in Cancer Subtyping

#' @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))
  }
}