R/KW.hit.express.R

In GRIN2: Genomic Random Interval (GRIN)

#' Associate Lesion Groups with Gene Expression
#'
#' @description
#' Performs the Kruskal-Wallis test to evaluate associations between genomic
#' lesion groups and expression levels of the corresponding genes.
#'
#' @usage
#' KW.hit.express(
#'   alex.data,
#'   gene.annotation,
#'   min.grp.size = NULL
#' )
#'
#' @param alex.data Output from the \code{\link{alex.prep.lsn.expr}} function.
#' A list containing:
#' \itemize{
#'   \item \code{alex.expr}: Gene expression data with genes represented by
#'   Ensembl gene IDs in rows and subjects in columns.
#'   \item \code{alex.lsn}: Lesion-group data for the same genes and subjects
#'   and in the same order as \code{alex.expr}.
#'   \item \code{alex.row.mtch}: A data frame containing the matched Ensembl
#'   gene IDs from the expression and lesion data.
#' }
#'
#' @param gene.annotation A gene annotation data frame, either supplied by the
#' user or retrieved using \code{\link{get.ensembl.annotation}}. It must contain
#' the columns \code{"gene"} (Ensembl gene ID), \code{"chrom"} (chromosome),
#' \code{"loc.start"} (gene start position), and \code{"loc.end"} (gene end
#' position).
#'
#' @param min.grp.size Optional numeric value specifying the minimum number of
#' subjects required in a lesion group for that group to be included in the
#' Kruskal-Wallis test. For a gene to be tested, at least two groups must each
#' contain at least \code{min.grp.size} subjects. Genes that do not satisfy this
#' requirement are retained in the output with an \code{NA} Kruskal-Wallis
#' p value.
#'
#' @details
#' For each matched gene, the function compares expression values across the
#' lesion groups defined in \code{alex.lsn} using the Kruskal-Wallis test.
#' Expression values are obtained from the corresponding gene in
#' \code{alex.expr}.
#'
#' Subjects without a genomic lesion affecting the gene are represented by the
#' lesion group \code{"none"}. Subjects affected by more than one lesion type
#' in the same gene may be represented by the group \code{"multiple"}, as
#' defined during preparation of the lesion matrix.
#'
#' In addition to the Kruskal-Wallis p value, the function reports the number
#' of subjects and the mean, median, and standard deviation of expression for
#' each lesion group represented in the input data.
#'
#' Kruskal-Wallis p values are adjusted for multiple testing using the
#' Benjamini-Hochberg false discovery rate procedure implemented by
#' \code{\link[stats]{p.adjust}}.
#'
#' @return
#' A data frame containing gene annotation information and lesion-expression
#' association results. Each row corresponds to a matched gene. Results include:
#' \itemize{
#'   \item \code{p.KW}: Kruskal-Wallis test p value.
#'   \item \code{q.KW}: FDR-adjusted q value.
#'   \item Columns ending in \code{_n.subjects}: number of subjects in each
#'   lesion group.
#'   \item Columns ending in \code{_mean}: mean expression within each lesion
#'   group.
#'   \item Columns ending in \code{_median}: median expression within each
#'   lesion group.
#'   \item Columns ending in \code{_sd}: standard deviation of expression within
#'   each lesion group.
#' }
#'
#' Genes that do not meet the minimum group-size requirement are retained with
#' \code{NA} values for \code{p.KW} and \code{q.KW}.
#'
#' @export
#'
#' @importFrom stats kruskal.test p.adjust
#'
#' @references
#' Kruskal, W. H., & Wallis, W. A. (1952). Use of Ranks in One-Criterion
#' Variance Analysis. Journal of the American Statistical Association,
#' 47(260), 583-621.
#'
#' Cao, X., Elsayed, A. H., & Pounds, S. B. (2023). Statistical Methods Inspired
#' by Challenges in Pediatric Cancer Multi-omics.
#'
#' @author
#' Abdelrahman Elsayed \email{abdelrahman.elsayed@stjude.org} and
#' Stanley Pounds \email{stanley.pounds@stjude.org}
#'
#' @seealso
#' \code{\link{alex.prep.lsn.expr}}
#'
#' @examples
#' data(expr_data)
#' data(lesion_data)
#' data(hg38_gene_annotation)
#'
#' # Prepare matched lesion-expression data (genes x subjects)
#' alex.data <- alex.prep.lsn.expr(expr_data,
#'                                 lesion_data,
#'                                 hg38_gene_annotation,
#'                                 min.expr = 1,
#'                                 min.pts.lsn = 5)
#'
#' # Perform Kruskal-Wallis test between lesion groups and expression levels
#' alex.kw.results <- KW.hit.express(alex.data,
#'                                   hg38_gene_annotation,
#'                                   min.grp.size = 5)
KW.hit.express=function(alex.data,          # output of the alex.prep.lsn.expr function (list of three data tables "alex.expr" with expression data ready for KW test, "alex.lsn" with lesion data and row.mtch)
                        gene.annotation,    # gene annotation data. First column "gene" should has ensembl IDs
                        min.grp.size=NULL)  # minimum group size to perform the test (there should be at least two groups with number of patients > min.grp.size)

{
  # Validate input data

  if (!is.list(alex.data))
    stop("alex.data must be the output from alex.prep.lsn.expr().")

  required.alex.objects=c("alex.expr","alex.lsn","alex.row.mtch")
  if (!all(required.alex.objects %in% names(alex.data)))
    stop("alex.data must contain: alex.expr, alex.lsn, and alex.row.mtch.")

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

  required.annotation.cols=c("gene","chrom","loc.start","loc.end")
  if (!all(required.annotation.cols %in% colnames(gene.annotation)))
    stop("gene.annotation must contain the columns: ",paste(required.annotation.cols,collapse=", "), ".")

  if (anyDuplicated(gene.annotation$gene))
    stop("The 'gene' column of gene.annotation must contain unique Ensembl gene IDs.")

  if (!is.null(min.grp.size) && (!is.numeric(min.grp.size) || length(min.grp.size)!=1 || is.na(min.grp.size) || min.grp.size<0))
    stop("min.grp.size must be NULL or a single non-negative numeric value.")

  expr.matrix=as.matrix(alex.data$alex.expr)
  lsn.matrix=as.matrix(alex.data$alex.lsn)
  row.mtch=alex.data$alex.row.mtch

  if (!identical(dim(expr.matrix),dim(lsn.matrix)))
    stop("alex.expr and alex.lsn must have the same dimensions.")

  if (!identical(rownames(expr.matrix),rownames(lsn.matrix)))
    stop("alex.expr and alex.lsn must contain genes in the same order.")

  if (!identical(colnames(expr.matrix),colnames(lsn.matrix)))
    stop("alex.expr and alex.lsn must contain subjects in the same order.")

  message(paste0("Computing KW P-value: ",date()))
  p=apply(row.mtch,1,one.KW.pvalue,
          expr.mtx=expr.matrix,
          hit.grps=lsn.matrix,
          min.grp.size=min.grp.size)

  message(paste0("Preparing results table: ",date()))
  kw.res=cbind(row.mtch,p.KW=p)

  # To prepare and print KW test results

  colnames(kw.res)=c("expr.row", "gene", "p.KW")
  kw.res.annotated=merge(gene.annotation,kw.res,by="gene", all.y=TRUE, sort=FALSE)
  kw.res.annotated$expr.row=NULL

  # Compute FDR adjusted q values

  kw.res.annotated$q.KW=stats::p.adjust(kw.res.annotated$p.KW,method="fdr")

  # To compute number of subjects, mean, median and stdev by lesion type

  unique.grps=sort(unique(lsn.matrix))
  unique.grps=sort(unique(unique.grps))
  count.by.lesion <- sapply(unique.grps,function(x)rowSums(lsn.matrix==x))
  colnames(count.by.lesion) = paste(colnames(count.by.lesion),"n.subjects",sep="_")
  count.by.lesion=as.matrix(count.by.lesion)

  mean.by.lsn=row.stats.by.group(expr.matrix, lsn.matrix, mean)
  colnames(mean.by.lsn) = paste(colnames(mean.by.lsn),"mean",sep="_")
  mean.by.lsn=as.matrix(mean.by.lsn)

  median.by.lsn=row.stats.by.group(expr.matrix, lsn.matrix, median)
  colnames(median.by.lsn) = paste(colnames(median.by.lsn),"median",sep="_")
  median.by.lsn=as.matrix(median.by.lsn)

  sd.by.lsn=row.stats.by.group(expr.matrix, lsn.matrix, sd)
  colnames(sd.by.lsn) = paste(colnames(sd.by.lsn),"sd",sep="_")
  sd.by.lsn=as.matrix(sd.by.lsn)

  # To make sure that the lesion-group statistics are in the same gene order as the KW results

  gene.order=match(kw.res.annotated$gene,rownames(count.by.lesion))

  if (anyNA(gene.order))
    stop("Gene IDs in the KW results could not be matched to the lesion-group statistics.")

  count.by.lesion=count.by.lesion[gene.order,,drop=FALSE]
  mean.by.lsn=mean.by.lsn[gene.order,,drop=FALSE]
  median.by.lsn=median.by.lsn[gene.order,,drop=FALSE]
  sd.by.lsn=sd.by.lsn[gene.order,,drop=FALSE]

  kw.res.final=cbind(kw.res.annotated, count.by.lesion, mean.by.lsn, median.by.lsn, sd.by.lsn)
  res=kw.res.final

  return(res)

}