R/alex.pathway.R

In GRIN2: Genomic Random Interval (GRIN)

#' Visualize Lesion and Expression Data by Pathway
#'
#' @description
#' Visualizes lesion and expression data for genes associated with a selected
#' biological pathway. Subjects are ordered by hierarchical clustering based on
#' their lesion profiles across pathway genes, and pathway genes are ordered
#' based on similarities in their lesion profiles.
#'
#' The function generates two aligned panels showing lesion and expression data
#' and returns the corresponding ordered data.
#'
#' @usage
#' alex.pathway(
#'   alex.data,
#'   lsn.clrs = NULL,
#'   lsn.data,
#'   pathways,
#'   selected.pathway
#' )
#'
#' @param alex.data Output from the \code{\link{alex.prep.lsn.expr}} function.
#' A list containing \code{"alex.expr"} (expression data), \code{"alex.lsn"}
#' (lesion-group data), and \code{"alex.row.mtch"} (matched Ensembl gene IDs).
#' The expression and lesion data contain genes in rows and subjects in columns,
#' with matching row and column order.
#'
#' @param lsn.clrs Optional named vector specifying colors for lesion groups.
#' Names must correspond to lesion types represented in the data. If
#' \code{NULL}, colors for individual lesion types are assigned using
#' \code{\link{default.grin.colors}}, with additional colors assigned to
#' \code{"none"} and \code{"multiple"} groups.
#'
#' @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).
#'
#' @param pathways A data frame containing pathway annotations with the columns
#' \code{"gene.name"} (gene symbol), \code{"ensembl.id"} (Ensembl gene ID),
#' and \code{"pathway"} (pathway name).
#'
#' @param selected.pathway A character string specifying the pathway to
#' visualize. The value must match a pathway in the \code{"pathway"} column of
#' \code{pathways}.
#'
#' @details
#' Genes annotated to \code{selected.pathway} are identified from
#' \code{pathways}, and their lesion and expression data are extracted from
#' \code{alex.data}.
#'
#' Pairwise lesion-profile distances are calculated using \code{dist.lsn}.
#' Hierarchical clustering is then used to order both pathway genes and subjects.
#' Subjects are clustered using lesion profiles across the selected pathway
#' genes. The same subject and gene ordering is used in both visualization
#' panels.
#'
#' The upper panel displays lesion groups using lesion-specific colors. The lower
#' panel displays standardized gene expression values, with lower expression
#' represented in blue, expression near the gene mean represented in white, and
#' higher expression represented in red.
#'
#' Gene symbols are used as row labels when available. If a gene symbol is
#' unavailable, the Ensembl gene ID is used instead.
#'
#' @return
#' A list containing:
#' \item{ordered.path.data}{A data frame containing lesion and expression data
#' for the selected pathway genes. Columns represent subjects in the order
#' determined by hierarchical clustering. Rows contain lesion data followed by
#' expression data and are labeled with \code{"_lsn"} and \code{"_expr"},
#' respectively.}
#'
#' The function also generates a figure with two aligned panels showing lesion
#' and expression data for the selected pathway genes. Both panels use the same
#' hierarchical clustering-based ordering.
#'
#' @export
#'
#' @importFrom tibble rownames_to_column
#' @importFrom stats hclust as.dist sd
#' @importFrom grDevices rgb
#' @importFrom graphics plot par rect legend text
#'
#' @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{alex.prep.lsn.expr}},
#' \code{\link[stats]{hclust}}
#'
#' @examples
#' data(expr_data)
#' data(lesion_data)
#' data(hg38_gene_annotation)
#' data(pathways)
#'
#' # Prepare matched expression and lesion data
#' alex.data <- alex.prep.lsn.expr(expr_data,
#'                                 lesion_data,
#'                                 hg38_gene_annotation,
#'                                 min.expr = 5,
#'                                 min.pts.lsn = 5)
#'
#' # Visualize pathway-level lesion and expression data using the JAK pathway
#' alex.path <- alex.pathway(alex.data,
#'                           lsn.data = lesion_data,
#'                           pathways = pathways,
#'                           selected.pathway = "Jak_Pathway")
#'
#' # Access the ordered data matrix used in the plot
#' alex.path$ordered.path.data[1:6,1:6]
alex.pathway=function(alex.data,          # output of the alex.prep.lsn.expr function (list of three data table "alex.expr" with expression data ready for KW test, "alex.lsn" with overlapped gene lesion data and row.mtch)
                      lsn.clrs=NULL,      # Specified colors per lesion types (gene plots when gene name is specified). If not specified, colors will be automatically assigned using default.grin.colors function
                      lsn.data,           # lesion data in a GRIN compatible format
                      pathways,           # data.table with three columns "gene.name" with gene symbols, "ensembl.id" with gene ensembl ID and "pathway" that has the pathway name
                      selected.pathway)   # pathway of interest

{

  # 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(lsn.data))
    stop("lsn.data must be a data frame.")

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

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

  required.pathway.cols=c("gene.name","ensembl.id","pathway")
  if (!all(required.pathway.cols %in% colnames(pathways)))
    stop("pathways must contain the columns: ",paste(required.pathway.cols,collapse=", "), ".")

  if (!is.character(selected.pathway) || length(selected.pathway)!=1 || is.na(selected.pathway) || selected.pathway=="")
    stop("selected.pathway must be a single pathway name.")

  if (!selected.pathway %in% pathways$pathway)
    stop("selected.pathway was not found in the 'pathway' column of pathways.")

  expr=alex.data$alex.expr
  lsn=alex.data$alex.lsn
  selected.genes=pathways[pathways$pathway==selected.pathway,]
  pathway.ensembl=selected.genes$ensembl.id

  # extract lesion and expression data for genes assigned to the pathway of interest

  path.expr=expr[(rownames(expr)%in%pathway.ensembl),,drop=FALSE]
  path.lsn=lsn[(rownames(lsn)%in%pathway.ensembl),,drop=FALSE]
  path.expr=as.matrix(path.expr)
  path.lsn=as.matrix(path.lsn)

  if (nrow(path.lsn)<2)
    stop("At least two genes from selected.pathway must be present in alex.data for hierarchical clustering.")

  if (!identical(dim(path.expr),dim(path.lsn)))
    stop("Pathway expression and lesion data must have the same dimensions.")

  if (!identical(rownames(path.expr),rownames(path.lsn)))
    stop("Pathway expression and lesion data must contain genes in the same order.")

  if (!identical(colnames(path.expr),colnames(path.lsn)))
    stop("Pathway expression and lesion data must contain subjects in the same order.")

  # assign colors to lesion groups

  lesions=lsn.data
  unique.grps=sort(unique(lesions$lsn.type))
  if (is.null(lsn.clrs))
  {
    lsn.grps.clr=default.grin.colors(unique.grps)
    common.grps.clr=c(none="gray", multiple="violet")
    lsn.clrs=c(lsn.grps.clr, common.grps.clr)
  }

  path.grps=unique(as.vector(path.lsn))
  if (!all(path.grps %in% names(lsn.clrs)))
    stop("lsn.clrs must provide colors for all lesion groups represented in the selected pathway.")

  lsn.grps=names(lsn.clrs)
  clrs=as.character(lsn.clrs)
  path.lsn.clr=path.lsn

  # to replace lesion groups with colors in the lsn matrix

  path.lsn.clr[path.lsn.clr %in% lsn.grps] <- clrs[match(path.lsn.clr, lsn.grps, nomatch = 0)]
  path.genes=selected.genes$ensembl.id
  path.gene.names=selected.genes$gene.name
  path.gene.names[is.na(path.gene.names) | path.gene.names==""]=path.genes[is.na(path.gene.names) | path.gene.names==""]
  names(path.gene.names)=path.genes

  # to replace ensembl IDs with gene name

  rownames(path.expr)=path.gene.names[rownames(path.expr)]
  rownames(path.lsn)=path.gene.names[rownames(path.lsn)]

  # compute the distance between each two genes based on the lesion data using dist.lsn function

  dist.lsn.genes=dist.lsn(t(path.lsn))
  hcl.lsn.genes=stats::hclust(stats::as.dist(dist.lsn.genes),"complete")

  # compute distance between subjects based on lesions affecting pathway genes using dist.lsn function

  path.lsn.dist=dist.lsn(path.lsn)
  path.lsn.dist=as.matrix(path.lsn.dist)

  path.ls.hcl=stats::hclust(stats::as.dist(path.lsn.dist),method="ward.D2")
  sum.none=rowSums(path.lsn.clr=="gray")
  subj.ord=path.ls.hcl$order
  subj.labels=path.ls.hcl$labels

  #####################################

  # side-by-side heatmap

  opar<-graphics::par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
  on.exit(graphics::par(opar), add = TRUE)
  plot(c(0,+1.1*ncol(path.lsn.clr)),
       c(0,-(nrow(path.lsn.clr)+1+nrow(path.expr))),
       type="n",xlab="",ylab="",
       axes=F)

  for (i in 1:nrow(path.lsn.clr))
    graphics::rect(0:(ncol(path.lsn.clr)-1),-(i-1),
                   1:ncol(path.lsn.clr),-i,
                   col=path.lsn.clr[hcl.lsn.genes$order[i],subj.ord],
                   border=NA)
  graphics::text(ncol(path.lsn),-(1:nrow(path.lsn))+0.5,
                 rownames(path.expr)[hcl.lsn.genes$order],
                 pos=4,cex=0.75)

  # Add legend

  lsn.inc=names(lsn.clrs)%in%path.lsn
  graphics::legend("topright", inset=c(-0.25,0.01),
                   fill=unlist(lsn.clrs[lsn.inc]),
                   legend=names(lsn.clrs[lsn.inc]),
                   cex=0.72,border=NA,bty="n")

  for (i in 1:nrow(path.expr))
  {
    y=path.expr[hcl.lsn.genes$order[i],]
    z=(y-mean(y))/stats::sd(y)
    clr=grDevices::rgb((z>0),0,(z<0),alpha=sqrt(1-exp(-abs(z))))
    graphics::rect(0:(ncol(path.lsn.clr)-1),-nrow(path.lsn.clr)-3-(i-1),
                   1:ncol(path.lsn.clr),-nrow(path.lsn.clr)-3-i,
                   col=clr[subj.ord],
                   border=NA)
  }

  graphics::legend("bottomright", inset=c(-0.25,0.2),
                   legend=c("Z.expr<0",0,"Z.expr>0"),
                   fill = c("blue", "white", "red"),
                   cex=0.72,border=NA,bty="n")

  graphics::text(ncol(path.expr),-nrow(path.lsn.clr)-3-1:nrow(path.expr)+0.5,
                 rownames(path.expr)[hcl.lsn.genes$order],cex=0.75,pos=4)

  # Extract ordered lesion and expression data for the pathway genes

  pts.labels=as.data.frame(subj.labels)
  pts.labels<-tibble::rownames_to_column(pts.labels, "subj.ord")
  pts.order=as.data.frame(subj.ord)
  pts.order$index=1:nrow(pts.order)
  ordered.subj.final=merge(pts.labels,pts.order,by="subj.ord", all.y=TRUE)
  ordered.subj.final=ordered.subj.final[order(ordered.subj.final$index),]
  path.lsn.df=as.data.frame(path.lsn)
  ordered.lsn.data=path.lsn.df[ordered.subj.final$subj.labels]
  rownames(ordered.lsn.data) = paste(rownames(ordered.lsn.data),"_lsn")
  path.expr.df=as.data.frame(path.expr)
  ordered.expr.data=path.expr.df[ordered.subj.final$subj.labels]
  rownames(ordered.expr.data) = paste(rownames(ordered.expr.data),"_expr")
  ordered.path.data=rbind(ordered.lsn.data, ordered.expr.data)

  res=list(ordered.path.data=ordered.path.data)

  return(res)

}