R/grin.logRank.R

In GRIN2: Genomic Random Interval (GRIN)

#' Log-Rank Test for Associations Between Genomic Lesions and Survival Outcomes
#'
#' @description
#' Performs gene-level log-rank tests to evaluate associations between genomic
#' lesions and time-to-event outcomes. For each gene, subjects are grouped
#' according to their lesion status and survival distributions are compared
#' using the log-rank test.
#'
#' @usage
#' grin.logRank(lsn.mtx,
#'              clin.data,
#'              annotation.data,
#'              clinvars,
#'              min.grp.size = NULL)
#'
#' @param lsn.mtx A gene-by-subject lesion matrix, typically generated using
#' \code{\link{prep.lsn.type.matrix}}. Rows represent genes and columns
#' represent subjects. For each gene-subject combination, entries indicate
#' lesion status: \code{"none"} if the gene is not affected, a specific lesion
#' type (e.g., \code{"mutation"} or \code{"fusion"}) if the gene is affected
#' by one lesion type, or \code{"multiple"} if the gene is affected by two or
#' more distinct lesion types in the same subject.
#'
#' @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{lsn.mtx}.
#'
#' @param annotation.data A gene annotation data frame containing a column
#' named \code{gene} with unversioned Ensembl gene IDs matching the gene IDs
#' used as row names in \code{lsn.mtx}. Annotation information is merged with
#' the final association results using these gene IDs.
#'
#' @param clinvars A character vector specifying the time-to-event clinical
#' variables to analyze. Each variable must be stored in \code{clin.data} as
#' a \code{\link[survival]{Surv}} object created using
#' \code{survival::Surv()}.
#'
#' @param min.grp.size Optional numeric value specifying the minimum number of
#' subjects required in each lesion-status group for a gene to be analyzed.
#' Groups may include \code{"none"}, individual lesion types such as
#' \code{"mutation"} or \code{"fusion"}, and \code{"multiple"}. Groups
#' containing fewer than \code{min.grp.size} subjects are excluded from the
#' analysis for that gene. A gene is tested only if at least two groups remain
#' after filtering.
#'
#' @details
#' Subject identifiers in the lesion matrix and clinical data are matched and
#' reordered before analysis.
#'
#' For each time-to-event outcome specified in \code{clinvars}, the function
#' performs a log-rank test for each gene using
#' \code{\link[survival]{survdiff}} with \code{rho = 0}.
#'
#' 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))}.
#'
#' The output also reports the number of subjects with and without an event
#' within each lesion group.
#'
#' @return
#' A data frame containing gene annotation information together with:
#' \itemize{
#'   \item \code{logRank_<endpoint>_pval}: Gene-level log-rank test p value
#'   for each survival endpoint specified in \code{clinvars}.
#'   \item \code{logRank_<endpoint>_qval}: Multiple-testing-adjusted q value
#'   for each survival endpoint.
#'   \item Numbers of subjects with an event in each lesion group.
#'   \item Numbers of subjects without an event in each lesion group.
#' }
#'
#' @export
#'
#' @references
#' Mantel, N. (1966). Evaluation of survival data and two new rank order
#' statistics arising in its consideration. Cancer Chemotherapy Reports,
#' 50(3), 163-170.
#'
#' Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate:
#' A practical and powerful approach to multiple testing. Journal of the Royal
#' Statistical Society: Series B, 57(1), 289-300.
#'
#' 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}
#'
#' @seealso \code{\link{prep.lsn.type.matrix}},
#' \code{\link[survival]{Surv}},
#' \code{\link[survival]{survdiff}}
#'
#' @examples
#' # Load the example datasets
#' data(lesion_data)
#' data(clin_data)
#' data(hg38_gene_annotation)
#'
#' # Prepare gene-level lesion data
#' gene.lsn <- prep.gene.lsn.data(lesion_data,
#'                                hg38_gene_annotation)
#'
#' # Identify overlaps between genomic lesions and genes
#' gene.lsn.overlap <- find.gene.lsn.overlaps(gene.lsn)
#'
#' # Create the lesion-type matrix
#' gene.lsn.type.mtx <- prep.lsn.type.matrix(gene.lsn.overlap,
#'                                           min.ngrp = 5)
#'
#' # Create the event-free survival object
#' clin_data$EFS <- survival::Surv(clin_data$efs.time,
#'                                 clin_data$efs.censor)
#'
#' # Run gene-level log-rank tests
#' logRank.efs <- grin.logRank(lsn.mtx = gene.lsn.type.mtx,
#'                             clin.data = clin_data,
#'                             annotation.data = hg38_gene_annotation,
#'                             clinvars = "EFS",
#'                             min.grp.size = 3)
#'
grin.logRank=function(lsn.mtx,
                      clin.data,
                      annotation.data,
                      clinvars,
                      min.grp.size=NULL)
{
  # Validate input data

  if (!is.matrix(lsn.mtx) && !is.data.frame(lsn.mtx))
    stop("lsn.mtx must be a matrix or 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 (anyDuplicated(clin.data$ID))
    stop("Subject IDs in clin.data$ID must be unique.")

  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(min.grp.size))
  {
    if (!is.numeric(min.grp.size) ||
        length(min.grp.size)!=1 ||
        is.na(min.grp.size) ||
        min.grp.size<1)
      stop("min.grp.size must be NULL or a positive numeric value.")
  }

  # Match and order subjects between lesion and clinical data

  lsn.mtx=t(lsn.mtx)
  lsn.df=as.data.frame(lsn.mtx)

  if (is.null(rownames(lsn.df)))
    stop("Subject IDs must be provided as column names of lsn.mtx.")

  if (anyDuplicated(rownames(lsn.df)))
    stop("Subject IDs in lsn.mtx must be unique.")

  clin.ids=as.character(clin.data$ID)
  common.ids=intersect(rownames(lsn.df),clin.ids)

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

  lsn.df=lsn.df[rownames(lsn.df) %in% common.ids,,drop=FALSE]
  clin.data=clin.data[clin.ids %in% common.ids,,drop=FALSE]

  lsn.df=lsn.df[order(rownames(lsn.df)),,drop=FALSE]
  clin.data=clin.data[order(as.character(clin.data$ID)),,drop=FALSE]

  if (!all(rownames(lsn.df)==as.character(clin.data$ID)))
    stop("Gene-lesion matrix subject IDs must match patient IDs in the clinical data.")

  lesion.columns=colnames(lsn.df)
  merged.data=cbind(lsn.df,clin.data)
  final.results=NULL

  # Run log-rank analysis for each survival endpoint

  for (v in seq_along(clinvars))
  {
    var.name=clinvars[v]
    thisvar=merged.data[[var.name]]

    # Confirm that the endpoint is a right-censored Surv object
    if (!survival::is.Surv(thisvar))
    {
      warning("Skipping ",var.name,
              ": grin.logRank() requires a survival outcome created with ",
              "survival::Surv().",call.=FALSE)
      next
    }

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

    message(paste0("Running log-rank test for association with ",
                   var.name,": ",date()))

    # Extract survival time and event status and remove missing outcomes
    surv.data=merged.data
    surv.data$surv.time=thisvar[,1]
    surv.data$surv.censor=thisvar[,2]

    keep=!is.na(surv.data$surv.time) & !is.na(surv.data$surv.censor)
    surv.data=surv.data[keep,,drop=FALSE]

    if (nrow(surv.data)==0)
    {
      warning("No subjects with non-missing survival data were available for ",
              var.name,".",call.=FALSE)
      next
    }

    if (!all(unique(surv.data$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
    }

    # Filter lesion groups according to the minimum group-size requirement
    surv.lsn.clms=surv.data[,lesion.columns,drop=FALSE]
    drop.cols=integer(0)

    for (i in seq_len(ncol(surv.lsn.clms)))
    {
      if (!is.null(min.grp.size))
      {
        group.counts=table(surv.lsn.clms[[i]],useNA="no")
        small.groups=names(group.counts[group.counts<min.grp.size])

        if (length(small.groups)>0)
          surv.lsn.clms[[i]][surv.lsn.clms[[i]] %in% small.groups]=NA
      }

      remaining.groups=unique(stats::na.omit(surv.lsn.clms[[i]]))

      if (length(remaining.groups)<2)
        drop.cols=c(drop.cols,i)
    }

    if (length(drop.cols)>0)
      surv.lsn.clms=surv.lsn.clms[,-drop.cols,drop=FALSE]

    if (ncol(surv.lsn.clms)==0)
    {
      message(paste0("No genes passed the lesion-group filtering for ",
                     var.name,"; skipping."))
      next
    }

    # Count subjects with and without events within each lesion group
    surv.event=surv.lsn.clms[surv.data$surv.censor==1,,drop=FALSE]
    surv.noevent=surv.lsn.clms[surv.data$surv.censor==0,,drop=FALSE]

    lesion.levels=sort(unique(c(unlist(surv.event,use.names=FALSE),
                                unlist(surv.noevent,use.names=FALSE))))
    lesion.levels=lesion.levels[!is.na(lesion.levels)]

    count.event=sapply(lesion.levels,function(x)
      colSums(surv.event==x,na.rm=TRUE))

    count.noevent=sapply(lesion.levels,function(x)
      colSums(surv.noevent==x,na.rm=TRUE))

    if (length(lesion.levels)==1)
    {
      count.event=matrix(count.event,ncol=1)
      count.noevent=matrix(count.noevent,ncol=1)
    }

    rownames(count.event)=colnames(surv.lsn.clms)
    rownames(count.noevent)=colnames(surv.lsn.clms)

    colnames(count.event)=paste0(lesion.levels,"_n.subjects.with.event")
    colnames(count.noevent)=paste0(lesion.levels,"_n.subjects.without.event")

    all.count=cbind(as.data.frame(count.event),
                    as.data.frame(count.noevent))

    # Perform gene-level log-rank tests

    pvalue.surv=rep(NA_real_,ncol(surv.lsn.clms))
    names(pvalue.surv)=colnames(surv.lsn.clms)

    for (i in seq_len(ncol(surv.lsn.clms)))
    {
      model.data=data.frame(surv.time=surv.data$surv.time,
                            surv.censor=surv.data$surv.censor,
                            lsn.group=surv.lsn.clms[[i]])

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

      if (nrow(model.data)==0 ||
          length(unique(model.data$lsn.group))<2)
        next

      fit=tryCatch(
        survival::survdiff(
          survival::Surv(surv.time,surv.censor)~lsn.group,
          data=model.data,
          rho=0,
          na.action=stats::na.omit),
        error=function(e) NULL)

      if (!is.null(fit))
        pvalue.surv[i]=fit$pvalue
    }

    # Compute FDR-adjusted q-values

    valid.p=!is.na(pvalue.surv)
    q.surv=rep(NA_real_,length(pvalue.surv))

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

    # Assemble results for this survival endpoint

    results.surv.final=data.frame(
      gene=names(pvalue.surv),
      pvalue.surv,
      q.surv,
      stringsAsFactors=FALSE)

    colnames(results.surv.final)[2:3]=c(
      paste0("logRank_",var.name,"_pval"),
      paste0("logRank_",var.name,"_qval"))

    all.count$gene=rownames(all.count)
    rownames(all.count)=NULL

    thisres=merge(results.surv.final,all.count,
                  by="gene",all.x=TRUE,sort=FALSE)

    # Distinguish lesion-group counts when multiple endpoints are analyzed
    if (length(clinvars)>1)
    {
      count.columns=grep("_n.subjects.",names(thisres),fixed=TRUE)
      names(thisres)[count.columns]=paste0(
        var.name,".",names(thisres)[count.columns])
    }

    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 survival outcomes were available for analysis.")

  # Add gene annotation information to the final results

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

  return(res.final)
}