R/decompose.gs.group.R

In mengchen18/mogsa: Multiple omics data integrative clustering and gene set analysis

#' Data-wise or PC-wise decomposition of gene set scores for all observations.
#' 
#' Data-wise or PC-wise decomposition of gene set scores (GSS) across all
#' observations. The predefined group/cluster information should be given so
#' that the mean decomposed GSSs for each group are returned and plotted.
#' 
#' This function could be used when the number of observation is large and
#' there are cluster/group information is available. In this case, the means of
#' decomposed gene set scores over each group is calculated. The vertical bar
#' on the end of each bar indicates the 95\% confident interval of the means.
#' 
#' @param x An object of class \code{\link{mgsa-class}} or
#' \code{\link{moa.sup-class}}
#' @param gs The gene set want to exam.
#' @param group An vector or factor to indicate the group of observations, such
#' as clusters.  See examples.
#' @param decomp A charater string either "data" or "pc" to indicate how the
#' gene set scores should be decomposed (with respect to data or PC.
#' @param nf The number of axes/PCs to be calculated and plotted.
#' @param x.legend Used to control the position of legends.
#' @param y.legend Used to control the position of legends.
#' @param plot A logical indicates if a plot should be drawn.
#' @param main The main title of plot.
#' @param \dots Other arguments passed to \code{barplot}.
#' @return Return nothing or a matrix depends on how argument \code{plot} is
#' set.
#' @author Chen Meng
#' @seealso See Also \code{\link{decompose.gs.ind}}
#' @references TBA
#' @importFrom matrixStats rowSds
#' @importFrom gplots barplot2
#' @export
#' @examples
#' 
#'   # library(mogsa)
#'   # loading gene expression data and supplementary data
#'   data(NCI60_4array_supdata)
#'   data(NCI60_4arrays)
#' 
#'   # using a list of data.frame as input
#'   mgsa <- mogsa(x = NCI60_4arrays, sup=NCI60_4array_supdata, nf=9,
#'                 proc.row = "center_ssq1", w.data = "inertia", statis = TRUE)
#' 
#'   colcode <- as.factor(sapply(strsplit(colnames(NCI60_4arrays$agilent), split="\\."), "[", 1))
#'   decompose.gs.group(x = mgsa, gs = 2, group = colcode, decomp = "data", plot = TRUE)
#'   decompose.gs.group(x = mgsa, gs = 2, group = colcode, decomp = "pc", nf = 3, plot = TRUE)
#' 
decompose.gs.group <-
function(x, gs, group, decomp = "data", nf=2, x.legend="bottomleft", y.legend=NULL, plot=TRUE, 
  main=NULL, ...) {
  
  # function only used inside function
  sedata <- function(ob) {
    m <- sapply(scl, function(x) 
      rowSums(sapply(x, function(y) y[gs, ob])))
    se <- rowSds(t(m))/sqrt(nrow(m))
    return(se)
  }
  
  sepc <- function(ob) {
    m <- sapply(scl, function(x) 
      sapply(x, function(y) y[gs, ob]))
    l <- split(m, rep(1:nf, rep(length(ob), nf)))
    l <- lapply(l, function(x) x[x!= 0])
    se <- sapply(l, sd)/sqrt(length(l[[1]]))
    names(se) <- names(scl[[1]])
    return(se)
  }
  # done
  
  if (inherits(x, "mgsa"))
    x <- x@sup
  scl <- lapply(x@score.sep, function (x) x[1:nf])
  
  if (is.list(group))
    cls <- group else {
      if (is.null(names(group)))
        names(group) <- paste("V", 1:length(group), sep = "")
        cls <- split(1:length(group), group)
    }
    
  gsm <- lapply(cls, function(ob) 
    decompose.gs.ind(x = x, plot = FALSE, gs=gs, obs = ob, nf = nf))
  
  if (decomp == "data") {
   s <- sapply(gsm, colSums) # by data
   sed <- sapply(cls, sedata)
  } else if (decomp == "pc") {
    s <- sapply(gsm, rowSums) # by pc
    sed <- sapply(cls, sepc)
  } else
    stop("unknown setting of decomp, decomp should be either data or pc!")
  
  if (plot) {
    ci.l <- s-1.96*sed
    ci.u <- s+1.96*sed
    col <- gray.colors(nrow(s))
    u <- barplot2(s, beside = TRUE, plot.ci = TRUE, ci.l = ci.l, ci.u = ci.u, col=col, 
      ylab = paste(decomp, "-wise decomposed gene set scores", sep=""), main=main)
    legend(x = x.legend, y = y.legend, legend = rownames(s), col=col, pch=15)
  }
  return(invisible(list(decomp.mean=s, decomp.se=sed)))
}