R/grin.assoc.expr.outcome.R

In GRIN2: Genomic Random Interval (GRIN)

#' Association Between Gene Expression and Clinical Outcomes
#'
#' @description
#' Performs gene-level association analyses between gene expression and
#' clinical outcomes. Time-to-event outcomes stored as \code{Surv} objects
#' are analyzed using Cox proportional hazards regression, while binary
#' outcomes coded as 0 and 1 are analyzed using logistic regression.
#'
#' For each clinical outcome, a separate regression model is fitted for every
#' gene to evaluate the association between that gene's expression and the
#' outcome. Models can optionally be adjusted for one or more clinical
#' covariates.
#'
#' @usage
#' grin.assoc.expr.outcome(expr.mtx,
#'                         clin.data,
#'                         annotation.data,
#'                         clinvars,
#'                         covariate = NULL)
#'
#' @param expr.mtx A data frame containing gene expression data with genes in
#' rows and subjects in columns. The first column must be named \code{"gene"}
#' and contain unique, unversioned Ensembl gene IDs (e.g.,
#' \code{"ENSG00000148400"}). Ensembl version suffixes such as \code{".5"} in
#' \code{"ENSG00000148400.5"} are not supported and should be removed before
#' analysis. Gene symbols should be converted to Ensembl gene IDs before using
#' this function. All remaining columns must correspond to subjects and contain
#' numeric expression values.
#'
#' @param clin.data A data frame containing clinical information. The data
#' frame must contain a column named \code{ID} with subject identifiers that
#' correspond to the subject identifiers in \code{expr.mtx}.
#'
#' @param annotation.data A gene annotation data frame containing a column
#' named \code{gene} with unversioned Ensembl gene IDs matching those in the
#' \code{gene} column of \code{expr.mtx}. Annotation information is merged with
#' the association results using these gene IDs.
#'
#' @param clinvars A character vector specifying the clinical outcome
#' variables to analyze. Time-to-event outcomes must be stored in
#' \code{clin.data} as \code{\link[survival]{Surv}} objects created using
#' \code{survival::Surv()}. Binary outcomes must be numeric and coded as
#' 0 and 1.
#'
#' @param covariate Optional character vector specifying one or more clinical
#' covariates to include in the regression models. Covariates may be
#' categorical or numeric. If \code{NULL}, models are fitted without
#' covariate adjustment.
#'
#' @details
#' Subject identifiers in the expression data and clinical data are matched
#' and reordered before analysis.
#'
#' For time-to-event outcomes stored as \code{Surv} objects, a separate Cox
#' proportional hazards model is fitted for each gene using
#' \code{\link[survival]{coxph}}. The reported hazard ratio represents the
#' relative change in hazard associated with a one-unit increase in the
#' expression value used as input to the model.
#'
#' Cox models that generate convergence warnings are excluded from downstream
#' inference. The hazard ratio, confidence interval, p value, and q value for
#' these genes are returned as \code{NA}, and a warning reports the number of
#' models that did not converge. When categorical covariates are included,
#' sparse categories or categories with few or no outcome events may lead to
#' model convergence problems. Users should review the distribution of
#' covariates and outcome events when a large number of adjusted models fail
#' to converge.
#'
#' For binary outcomes coded numerically as 0 and 1, a separate logistic
#' regression model is fitted for each gene using \code{\link[stats]{glm}}
#' with \code{family = "binomial"}. The reported odds ratio represents the
#' change in the odds of the outcome coded as 1 associated with a one-unit
#' increase in expression.
#'
#' Numeric outcomes containing values other than 0 and 1 are not analyzed
#' and generate a warning. Continuous outcomes should be analyzed using an
#' appropriate regression model outside this function.
#'
#' When \code{covariate} is provided, the specified covariates are included
#' together with gene expression in each Cox proportional hazards or logistic
#' regression model.
#'
#' P values are adjusted for multiple testing using the Benjamini-Hochberg
#' false discovery rate procedure together with the Pounds and Cheng estimator
#' of the proportion of tests having a true null hypothesis:
#' \code{pi.hat = min(1, 2 * mean(p))}.
#'
#' Cox proportional hazards models assume proportional hazards. Users should
#' evaluate this assumption when interpreting genes of particular interest.
#'
#' @return
#' A data frame containing gene annotation information together with
#' outcome-specific association results.
#'
#' For time-to-event outcomes, results include:
#' \itemize{
#'   \item Hazard ratio (HR).
#'   \item Lower bound of the 95 percent confidence interval.
#'   \item Upper bound of the 95 percent confidence interval.
#'   \item Cox proportional hazards model p value.
#'   \item Multiple-testing-adjusted q value.
#' }
#'
#' For binary outcomes, results include:
#' \itemize{
#'   \item Odds ratio (OR).
#'   \item Lower bound of the 95 percent confidence interval.
#'   \item Upper bound of the 95 percent confidence interval.
#'   \item Logistic regression p value.
#'   \item Multiple-testing-adjusted q value.
#' }
#'
#' Cox models that do not converge are retained in the output with
#' \code{NA} values for their model statistics.
#'
#' @export
#'
#' @references
#' Cox, D. R. (1972). Regression Models and Life-Tables.
#' Journal of the Royal Statistical Society: Series B (Methodological),
#' 34(2), 187-202.
#'
#' Nelder, J. A., & Wedderburn, R. W. M. (1972). Generalized Linear Models.
#' Journal of the Royal Statistical Society: Series A (General),
#' 135(3), 370-384.
#'
#' Pounds, S., & Cheng, C. (2006). Robust estimation of the false discovery
#' rate. Bioinformatics, 22(16), 1979-1987.
#'
#' @author
#' Abdelrahman Elsayed \email{abdelrahman.elsayed@stjude.org} and
#' Stanley Pounds \email{stanley.pounds@stjude.org}
#'
#' @examples
#' # Load the example datasets
#' data(expr_data)
#' data(clin_data)
#' data(hg38_gene_annotation)
#'
#' # Create the event-free survival object
#' clin_data$EFS <- survival::Surv(clin_data$efs.time,
#'                                 clin_data$efs.censor)
#'
#' # Specify the survival endpoint
#' clinvars <- c("EFS")
#'
#' # Run Cox proportional hazards models
#' coxph.efs <- grin.assoc.expr.outcome(
#'   expr.mtx = expr_data,
#'   clin.data = clin_data,
#'   annotation.data = hg38_gene_annotation,
#'   clinvars = clinvars
#' )
#'
#' # Run Cox proportional hazards models with covariate adjustment
#' coxph.efs.adj <- grin.assoc.expr.outcome(
#'   expr.mtx = expr_data,
#'   clin.data = clin_data,
#'   annotation.data = hg38_gene_annotation,
#'   clinvars = clinvars,
#'   covariate = "WBC"
#' )
#'
#' # A binary outcome coded as 0 and 1 can also be analyzed
#' clinvars <- c("MRD.binary")
#'
#' logistic.mrd <- grin.assoc.expr.outcome(
#'   expr.mtx = expr_data,
#'   clin.data = clin_data,
#'   annotation.data = hg38_gene_annotation,
#'   clinvars = clinvars
#' )
#'
grin.assoc.expr.outcome <- function(expr.mtx,
                                    clin.data,
                                    annotation.data,
                                    clinvars,
                                    covariate = NULL)
{
  # Validate input data
  if (!is.data.frame(expr.mtx))
    stop("expr.mtx must be a data frame.")

  if (ncol(expr.mtx) < 2)
    stop("expr.mtx must contain a 'gene' column and at least one subject expression column.")

  if (names(expr.mtx)[1] != "gene")
    stop("The first column of expr.mtx must be named 'gene'.")

  if (!is.data.frame(clin.data))
    stop("clin.data must be a data frame.")

  if (!is.data.frame(annotation.data))
    stop("annotation.data must be a data frame.")

  if (!"ID" %in% names(clin.data))
    stop("clin.data must contain a column named 'ID'.")

  if (!"gene" %in% names(annotation.data))
    stop("annotation.data must contain a column named 'gene'.")

  if (length(clinvars) == 0)
    stop("At least one clinical outcome must be specified in clinvars.")

  if (!all(clinvars %in% names(clin.data)))
  {
    missing.vars <- clinvars[!clinvars %in% names(clin.data)]
    stop("The following clinvars were not found in clin.data: ",
         paste(missing.vars, collapse = ", "))
  }

  if (!is.null(covariate))
  {
    if (!all(covariate %in% names(clin.data)))
    {
      missing.covariates <- covariate[!covariate %in% names(clin.data)]
      stop("The following covariates were not found in clin.data: ",
           paste(missing.covariates, collapse = ", "))
    }

    if (any(covariate %in% clinvars))
      stop("Variables specified in clinvars cannot also be used as covariates.")
  }

  # Prepare and validate expression data

  gene.names <- as.character(expr.mtx$gene)

  if (anyNA(gene.names) || any(gene.names == ""))
    stop("The 'gene' column of expr.mtx must not contain missing or empty gene IDs.")

  if (anyDuplicated(gene.names))
    stop("The 'gene' column of expr.mtx must contain unique Ensembl gene IDs.")

  if (any(grepl("\\.[0-9]+$", gene.names)))
    stop("Ensembl gene IDs in expr.mtx must not include version suffixes (e.g., use 'ENSG00000148400' instead of 'ENSG00000148400.5'). Please remove Ensembl version suffixes before running the analysis.")

  if (!all(grepl("^ENSG[0-9]+$", gene.names)))
    stop("The 'gene' column of expr.mtx must contain Ensembl gene IDs (e.g., 'ENSG00000148400'). Gene symbols are not supported and should be converted to Ensembl gene IDs before analysis.")

  expr.mtx <- expr.mtx[, -1, drop = FALSE]

  if (anyDuplicated(colnames(expr.mtx)))
    stop("Subject IDs in expr.mtx must be unique.")

  if (anyDuplicated(clin.data$ID))
    stop("Subject IDs in clin.data$ID must be unique.")

  # Convert expression values to numeric

  original.na <- is.na(expr.mtx)
  expr.numeric <- suppressWarnings(
    matrix(as.numeric(as.matrix(expr.mtx)),
           nrow = nrow(expr.mtx), ncol = ncol(expr.mtx),
           dimnames = list(gene.names, colnames(expr.mtx)))
  )

  invalid.values <- is.na(expr.numeric) & !original.na

  if (any(invalid.values))
    stop("Expression values could not be converted to numeric. ",
         "Please check expr.mtx for non-numeric values.")

  expr.mtx <- expr.numeric

  # Match subjects between expression and clinical data

  clin.ids <- as.character(clin.data$ID)
  common.ids <- intersect(colnames(expr.mtx), clin.ids)

  if (length(common.ids) == 0)
    stop("No matching subject IDs were found between expr.mtx and clin.data.")

  expr.mtx <- expr.mtx[, common.ids, drop = FALSE]
  clin.data <- clin.data[match(common.ids, as.character(clin.data$ID)), , drop = FALSE]

  if (!all(colnames(expr.mtx) == as.character(clin.data$ID)))
    stop("Expression matrix subject IDs must match patient IDs ",
         "in the clinical data.")

  final.results <- NULL

  # Analyze each clinical outcome

  for (i in seq_along(clinvars))
  {
    var.name <- clinvars[i]
    thisvar <- clin.data[[var.name]]

    # Cox proportional hazards regression for survival outcomes

    if (survival::is.Surv(thisvar))
    {
      if (ncol(thisvar) != 2)
      {
        warning("Skipping ", var.name,
                ": grin.assoc.expr.outcome() currently supports right-censored ",
                "Surv(time, event) outcomes.", call. = FALSE)
        next
      }

      message(paste0("Running Cox proportional hazards models for association with ",
                     var.name, ": ", date()))

      surv.time <- thisvar[, 1]
      surv.censor <- thisvar[, 2]

      if (!all(unique(stats::na.omit(surv.censor)) %in% c(0, 1)))
      {
        warning("Skipping ", var.name,
                ": the event indicator must use 0 for censoring and 1 for events.",
                call. = FALSE)
        next
      }

      # Initialize Cox model results

      model.results <- matrix(
        NA_real_, nrow = nrow(expr.mtx), ncol = 4,
        dimnames = list(gene.names, c("HR", "lower95", "upper95", "pvalue"))
      )

      n.nonconverged <- 0

      # Fit one Cox model per gene

      for (g in seq_len(nrow(expr.mtx)))
      {
        model.data <- data.frame(
          time = surv.time, event = surv.censor,
          expression = as.numeric(expr.mtx[g, ]),
          stringsAsFactors = FALSE
        )

        if (!is.null(covariate))
          model.data <- cbind(model.data, clin.data[, covariate, drop = FALSE])

        model.data <- model.data[stats::complete.cases(model.data), , drop = FALSE]

        if (nrow(model.data) < 2)
          next

        if (length(unique(model.data$expression)) < 2)
          next

        if (sum(model.data$event == 1) < 1)
          next

        # Construct unadjusted or covariate-adjusted model formula

        if (is.null(covariate))
        {
          model.formula <- survival::Surv(time, event) ~ expression
        } else {
          model.formula <- stats::reformulate(
            c("expression", covariate),
            response = "survival::Surv(time, event)"
          )
        }

        # Fit Cox model and identify convergence warnings

        fit.warning <- FALSE
        fit <- tryCatch(
          withCallingHandlers(
            survival::coxph(model.formula, data = model.data),
            warning = function(w)
            {
              warning.message <- conditionMessage(w)

              if (grepl(paste("did not converge",
                              "ran out of iterations",
                              "coefficient may be infinite",
                              sep = "|"),
                        warning.message, ignore.case = TRUE))
              {
                fit.warning <<- TRUE
                invokeRestart("muffleWarning")
              }
            }
          ),
          error = function(e) NULL
        )

        if (fit.warning)
        {
          n.nonconverged <- n.nonconverged + 1
          next
        }

        if (is.null(fit))
          next

        fit.summary <- summary(fit)

        if (!"expression" %in% rownames(fit.summary$coefficients))
          next

        model.results[g, "HR"] <-
          fit.summary$coefficients["expression", "exp(coef)"]

        model.results[g, "lower95"] <-
          fit.summary$conf.int["expression", "lower .95"]

        model.results[g, "upper95"] <-
          fit.summary$conf.int["expression", "upper .95"]

        model.results[g, "pvalue"] <-
          fit.summary$coefficients["expression", "Pr(>|z|)"]
      }

      if (n.nonconverged > 0)
        warning(n.nonconverged, " Cox proportional hazards model(s) for ",
                var.name, " did not converge and were returned as NA.",
                call. = FALSE)

      # Compute FDR-adjusted q-values

      pvalue <- model.results[, "pvalue"]
      valid.p <- !is.na(pvalue)
      qvalue <- rep(NA_real_, length(pvalue))

      if (any(valid.p))
      {
        pi.hat <- min(1, 2 * mean(pvalue[valid.p]))
        qvalue[valid.p] <- pi.hat *
          stats::p.adjust(pvalue[valid.p], method = "fdr")
        qvalue[qvalue > 1] <- 1
      }

      # Assemble Cox results

      thisres <- data.frame(
        gene = gene.names,
        HR = model.results[, "HR"],
        lower95 = model.results[, "lower95"],
        upper95 = model.results[, "upper95"],
        pvalue = pvalue,
        qvalue = qvalue,
        stringsAsFactors = FALSE
      )

      if (is.null(covariate))
      {
        colnames(thisres)[2:6] <- c(
          paste0("cox_", var.name, "_HR"),
          paste0("cox_", var.name, "_lower95"),
          paste0("cox_", var.name, "_upper95"),
          paste0("cox_", var.name, "_pval"),
          paste0("cox_", var.name, "_qval")
        )
      } else {
        colnames(thisres)[2:6] <- c(
          paste0("cox_", var.name, "_HR_adj"),
          paste0("cox_", var.name, "_lower95_adj"),
          paste0("cox_", var.name, "_upper95_adj"),
          paste0("cox_", var.name, "_pval_adj"),
          paste0("cox_", var.name, "_qval_adj")
        )
      }

      # Logistic regression for binary outcomes

    } else if (is.numeric(thisvar) &&
               all(stats::na.omit(unique(thisvar)) %in% c(0, 1)))
    {
      message(paste0("Running logistic regression models for association with ",
                     var.name, ": ", date()))

      # Initialize logistic model results

      model.results <- matrix(
        NA_real_, nrow = nrow(expr.mtx), ncol = 4,
        dimnames = list(gene.names, c("OR", "lower95", "upper95", "pvalue"))
      )

      # Fit one logistic regression model per gene

      for (g in seq_len(nrow(expr.mtx)))
      {
        model.data <- data.frame(
          outcome = thisvar,
          expression = as.numeric(expr.mtx[g, ]),
          stringsAsFactors = FALSE
        )

        if (!is.null(covariate))
          model.data <- cbind(model.data, clin.data[, covariate, drop = FALSE])

        model.data <- model.data[stats::complete.cases(model.data), , drop = FALSE]

        if (nrow(model.data) < 2)
          next

        if (length(unique(model.data$expression)) < 2)
          next

        if (length(unique(model.data$outcome)) < 2)
          next

        # Construct unadjusted or covariate-adjusted model formula

        if (is.null(covariate))
        {
          model.formula <- outcome ~ expression
        } else {
          model.formula <- stats::reformulate(
            c("expression", covariate),
            response = "outcome"
          )
        }

        fit <- tryCatch(
          stats::glm(model.formula, data = model.data,
                     family = stats::binomial()),
          error = function(e) NULL
        )

        if (is.null(fit))
          next

        fit.summary <- summary(fit)

        if (!"expression" %in% rownames(fit.summary$coefficients))
          next

        coefficient <- fit.summary$coefficients["expression", "Estimate"]
        standard.error <- fit.summary$coefficients["expression", "Std. Error"]

        model.results[g, "OR"] <- exp(coefficient)
        model.results[g, "lower95"] <- exp(coefficient - 1.96 * standard.error)
        model.results[g, "upper95"] <- exp(coefficient + 1.96 * standard.error)
        model.results[g, "pvalue"] <-
          fit.summary$coefficients["expression", "Pr(>|z|)"]
      }

      # Compute FDR-adjusted q-values

      pvalue <- model.results[, "pvalue"]
      valid.p <- !is.na(pvalue)
      qvalue <- rep(NA_real_, length(pvalue))

      if (any(valid.p))
      {
        pi.hat <- min(1, 2 * mean(pvalue[valid.p]))
        qvalue[valid.p] <- pi.hat *
          stats::p.adjust(pvalue[valid.p], method = "fdr")
        qvalue[qvalue > 1] <- 1
      }

      # Assemble logistic regression results

      thisres <- data.frame(
        gene = gene.names,
        odds.ratio = model.results[, "OR"],
        lower95 = model.results[, "lower95"],
        upper95 = model.results[, "upper95"],
        pvalue = pvalue,
        qvalue = qvalue,
        stringsAsFactors = FALSE
      )

      if (is.null(covariate))
      {
        colnames(thisres)[2:6] <- c(
          paste0("logistic_", var.name, "_OR"),
          paste0("logistic_", var.name, "_lower95"),
          paste0("logistic_", var.name, "_upper95"),
          paste0("logistic_", var.name, "_pval"),
          paste0("logistic_", var.name, "_qval")
        )
      } else {
        colnames(thisres)[2:6] <- c(
          paste0("logistic_", var.name, "_OR.adj"),
          paste0("logistic_", var.name, "_lower95.adj"),
          paste0("logistic_", var.name, "_upper95.adj"),
          paste0("logistic_", var.name, "_pval.adj"),
          paste0("logistic_", var.name, "_qval.adj")
        )
      }

    } else if (is.numeric(thisvar))
    {
      warning("Skipping ", var.name,
              ": numeric outcomes must be binary and coded as 0/1 ",
              "for logistic regression.", call. = FALSE)
      next

    } else {
      warning("Skipping ", var.name,
              ": outcomes must either be a survival object created with ",
              "survival::Surv() or a numeric binary variable coded as 0/1.",
              call. = FALSE)
      next
    }

    # Merge results across clinical outcomes

    if (is.null(final.results))
    {
      final.results <- thisres
    } else {
      final.results <- merge(final.results, thisres,
                             by = "gene", all = TRUE, sort = FALSE)
    }
  }

  if (is.null(final.results))
    stop("No valid clinical outcomes were available for analysis.")

  # Add gene annotation information

  res.final <- merge(annotation.data, final.results,
                     by = "gene", all.y = TRUE, sort = FALSE)

  return(res.final)
}