R/alex.prep.lsn.expr.R

In GRIN2: Genomic Random Interval (GRIN)

#' Prepare Lesion and Expression Data for Kruskal-Wallis Test
#'
#' @description
#' Prepares matched gene-level lesion and expression data matrices for use with
#' \code{\link{KW.hit.express}}, which performs Kruskal-Wallis tests to assess
#' associations between genomic lesion groups and gene expression levels.
#'
#' @usage
#' alex.prep.lsn.expr(
#'   expr.mtx,
#'   lsn.data,
#'   gene.annotation,
#'   min.expr = NULL,
#'   min.pts.lsn = NULL
#' )
#'
#' @param expr.mtx A data frame containing normalized 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 measurements; missing values (\code{NA}) are allowed.
#'
#' @param lsn.data A data frame containing genomic lesion data in GRIN-compatible
#' format. It must contain the columns \code{"ID"} (subject ID), \code{"chrom"}
#' (chromosome), \code{"loc.start"} (lesion start position), \code{"loc.end"}
#' (lesion end position), and \code{"lsn.type"} (lesion type; e.g., mutation,
#' gain, deletion, or structural rearrangement).
#'
#' @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.expr Optional numeric value specifying the minimum total expression
#' required for a gene to be retained. Total expression is calculated as the sum
#' of expression values across the matched subjects. If \code{NULL}, genes are
#' not filtered according to total expression.
#'
#' @param min.pts.lsn Optional numeric value specifying the minimum number of
#' subjects with at least one genomic lesion affecting a gene for that gene to
#' be retained. If \code{NULL}, genes are not filtered according to lesion
#' frequency.
#'
#' @details
#' Genomic lesions are first mapped to genes using
#' \code{\link{prep.gene.lsn.data}} and \code{\link{find.gene.lsn.overlaps}}.
#' The resulting gene-lesion overlaps are then converted to a gene-by-subject
#' lesion-type matrix using \code{\link{prep.lsn.type.matrix}}.
#'
#' Only subjects represented in both the lesion and expression datasets are
#' retained. Genes may additionally be filtered according to \code{min.expr}
#' and \code{min.pts.lsn}. The lesion and expression matrices are subsequently
#' restricted to their shared set of genes and ordered identically by Ensembl
#' gene ID and subject ID.
#'
#' The resulting matched matrices can be supplied directly to
#' \code{\link{KW.hit.express}} for lesion-expression association analysis.
#'
#' @return
#' A list with the following components:
#' \item{alex.expr}{A numeric matrix of gene expression data with Ensembl gene
#' IDs as row names and subject IDs as column names.}
#' \item{alex.lsn}{A data frame containing lesion groups for the same genes and
#' subjects and in the same order as \code{alex.expr}.}
#' \item{alex.row.mtch}{A data frame containing the matched Ensembl gene IDs
#' from the expression and lesion matrices.}
#'
#' @export
#'
#' @references
#' 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},
#' Stanley Pounds \email{stanley.pounds@stjude.org}
#'
#' @seealso
#' \code{\link{KW.hit.express}},
#' \code{\link{prep.gene.lsn.data}},
#' \code{\link{find.gene.lsn.overlaps}},
#' \code{\link{prep.lsn.type.matrix}},
#' \code{\link{get.ensembl.annotation}}
#'
#' @examples
#' data(expr_data)
#' data(lesion_data)
#' data(hg38_gene_annotation)
#'
#' # Prepare matched lesion and expression data
#' alex.data <- alex.prep.lsn.expr(expr_data,
#'                                 lesion_data,
#'                                 hg38_gene_annotation,
#'                                 min.expr = 1,
#'                                 min.pts.lsn = 5)
alex.prep.lsn.expr=function(expr.mtx,          # Normalized log2 transformed expression data with genes in rows and subjects in columns (first column "gene" should contain unversioned Ensembl gene IDs)
                            lsn.data,          # lesion data in GRIN compatible format. Data frame should has five columns that include "ID" with patient ID, "chrom" which is the chromosome on which the lesion is located, "loc.start" which is the lesion start position, "loc.end" the lesion end position and "lsn.type" which is the lesion type for example gain, loss, mutation, fusion, etc...
                            gene.annotation,   # gene annotation data either provided by the user or directly retreived using get.ensembl.annotation function included in the GRIN2 library. Object should has four columns "gene" which is the ensembl ID of annotated genes, "chrom" which is the chromosome on which the gene is located, "loc.start" which is the gene start position, and "loc.end" the gene end position.
                            min.expr=NULL,     # minimum allowed expression level of the gene (the sum of expression level of the gene in all patients; helpful to exclude genes with very low expression)
                            min.pts.lsn=NULL)  # minimum number of patients with any type of lesions in a certain gene

{
  # 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 (colnames(expr.mtx)[1]!="gene")
    stop("The first column of expr.mtx must be named 'gene'.")

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

  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 (!is.null(min.expr) && (!is.numeric(min.expr) || length(min.expr)!=1 || is.na(min.expr)))
    stop("min.expr must be NULL or a single numeric value.")

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

  # Prepare expression data for ALEX analysis

  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.")

  original.na=is.na(expr.mtx)
  expr.numeric=as.data.frame(lapply(expr.mtx,function(x) suppressWarnings(as.numeric(as.character(x)))),check.names=FALSE)
  rownames(expr.numeric)=gene.names
  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

  gene.lsn=prep.gene.lsn.data(lsn.data, gene.annotation)    # Prepare gene and lesion data for later computations
  gene.lsn.overlap= find.gene.lsn.overlaps(gene.lsn)        # Use the results of prep.gene.lsn.data to find lesion-gene overlaps

  message(paste0("Preparing gene-lesion matrix: ",date()))
  lsn.type.matrix=prep.lsn.type.matrix(gene.lsn.overlap)    # prepare lesion type matrix (just one row for each gene with all lesion types included)
  lsn.grp.mtx=as.data.frame(lsn.type.matrix)
  id.hits = colnames(lsn.grp.mtx)

  id.expr=colnames(expr.mtx)
  shared.ids=intersect(colnames(lsn.grp.mtx),id.expr)

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

  lsn.grp.mtx=lsn.grp.mtx[,shared.ids,drop=FALSE]
  expr.mtx=expr.mtx[,shared.ids,drop=FALSE]

  if (is.numeric(min.expr))
  {
    # To keep only genes with total expression level > value specified by the user in the min.expr argument
    total.expr=rowSums(expr.mtx,na.rm=TRUE)
    expr.grp.keep=total.expr>=min.expr
    expr.mtx=expr.mtx[expr.grp.keep,,drop=FALSE]
  }
  ############################

  # To exclude any gene that has no lesions in at least the number of patients specified in min.pts.lsn

  if (is.numeric(min.pts.lsn))
  {
    n.none=rowSums(lsn.grp.mtx=="none",na.rm=TRUE)
    lsn.grp.keep=(ncol(lsn.grp.mtx)-n.none)>=min.pts.lsn
    lsn.grp.mtx=lsn.grp.mtx[lsn.grp.keep,,drop=FALSE]
  }

  # To keep only shared set of genes between expression and lesion matrices

  lsn.genes=rownames(lsn.grp.mtx)
  expr.mtx=expr.mtx[rownames(expr.mtx) %in% lsn.genes,,drop=FALSE]

  expr.genes=rownames(expr.mtx)
  lsn.grp.mtx=lsn.grp.mtx[rownames(lsn.grp.mtx) %in% expr.genes,,drop=FALSE]

  if (nrow(expr.mtx)==0 || nrow(lsn.grp.mtx)==0)
    stop("No shared genes remain between the expression and lesion matrices after filtering.")

  ## order lsn matrix by gene ID first then colnames for patient IDs alphabetically

  lsn.grp.mtx=lsn.grp.mtx[order(rownames(lsn.grp.mtx)),,drop=FALSE]
  lsn.grp.mtx=lsn.grp.mtx[,order(colnames(lsn.grp.mtx)),drop=FALSE]

  ## order expr matrix by gene ID first then colnames for patient IDs alphabetically

  expr.mtx=expr.mtx[order(rownames(expr.mtx)),,drop=FALSE]
  expr.mtx=expr.mtx[,order(colnames(expr.mtx)),drop=FALSE]

  # To make sure that both lesion and expression matrices have same patients and genes order

  check.rownames=all(rownames(lsn.grp.mtx)==rownames(expr.mtx))
  if (check.rownames==FALSE)
    stop("Gene-lesion matrix rownames must match expression matrix rownames.")

  check.colnames=all(colnames(lsn.grp.mtx)==colnames(expr.mtx))
  if (check.colnames==FALSE)
    stop("Gene-lesion matrix patient IDs must match expression matrix patient IDs.")

  # To generate row.mtch file

  lsn.ensembl.id=rownames(lsn.grp.mtx)
  expr.ensembl.id=rownames(expr.mtx)

  row.mtch=cbind.data.frame(expr.row=expr.ensembl.id,
                            hit.row=lsn.ensembl.id)

  res=list(alex.expr=expr.mtx,      # expression data (data table with gene ensembl IDs as row names and each column is a patient)
           alex.lsn=lsn.grp.mtx,    # overlapped gene lesion data table with gene ensembl IDs as row names and each column is a patient
           alex.row.mtch=row.mtch)  # ensembl ID for one row of the alex.lsn table to be paired with one row of the expression data for the Kruskal-Wallis test

  return(res)

}