R/BayesSpikeSolver.R

In TReNA: Fit transcriptional regulatory networks using gene expression, priors, machine learning

#----------------------------------------------------------------------------------------------------
#' An S4 class to represent a Bayes Spike solver
#'
#' @import vbsr
#' @import methods
#' 
#' @include Solver.R
#' @name BayesSpikeSolver-class

.BayesSpikeSolver <- setClass ("BayesSpikeSolver", contains="Solver")
#----------------------------------------------------------------------------------------------------
#' Create a Solver class object using the Bayes Spike Solver
#' 
#' @param mtx.assay An assay matrix of gene expression data
#' @param quiet A logical denoting whether or not the solver should print output
#'
#' @export
#' 
#' @return A Solver class object with Bayes Spike as the solver
#'
#' @family Solver class objects
#' 
#' @seealso  \code{\link{solve.BayesSpike}}, \code{\link{getAssayData}}
#'
#' @examples
#' solver <- BayesSpikeSolver()

BayesSpikeSolver <- function(mtx.assay=matrix(), quiet=TRUE)
{
    obj <- .BayesSpikeSolver(Solver(mtx.assay=mtx.assay, quiet=quiet))

    # Send a warning if there's a row of zeros
    if(!is.na(max(mtx.assay)) & any(rowSums(mtx.assay) == 0))
       warning("One or more gene has zero expression; this may cause difficulty when using Bayes Spike. You may want to try 'lasso' or 'ridge' instead.")

    obj

} # BayesSpikeSolver, the constructor
#----------------------------------------------------------------------------------------------------
#' Run the Bayes Spike Solver
#'
#' @rdname solve.BayesSpike
#' @aliases run.BayesSpikeSolver solve.BayesSpike
#' 
#' @description Given a TReNA object with Bayes Spike as the solver, use the \code{\link{vbsr}}
#' function to estimate coefficients for each transcription factor as a predictor of the target
#' gene's expression level. This method should be called using the \code{\link{solve}} method on an
#' appropriate TReNA object. 
#'
#' @param obj An object of the class Solver with "bayesSpike" as the solver string
#' @param target.gene A designated target gene that should be part of the mtx.assay data
#' @param tfs The designated set of transcription factors that could be associated with the target gene.
#' @param tf.weights A set of weights on the transcription factors (default = rep(1, length(tfs)))
#' @param extraArgs Modifiers to the Bayes Spike solver; this includes \code{n_orderings}, the
#' number of random starts used by the solver
#' 
#' @return A data frame containing the coefficients relating the target gene to each transcription factor,
#' plus other fit parameters
#'
#' @seealso \code{\link{vbsr}}, \code{\link{BayesSpikeSolver}}
#'
#' @family solver methods
#' 
#' @examples
#' # Load included Alzheimer's data, create a TReNA object with Bayes Spike as solver, and solve
#' load(system.file(package="TReNA", "extdata/ampAD.154genes.mef2cTFs.278samples.RData"))
#' trena <- TReNA(mtx.assay = mtx.sub, solver = "bayesSpike")
#' target.gene <- "MEF2C"
#' tfs <- setdiff(rownames(mtx.sub), target.gene)
#' tbl <- solve(trena, target.gene, tfs)
#'
#' # Solve the same Alzheimer's problem, but this time set the number of random starts to 100
#' tbl <- solve(trena, target.gene, tfs, extraArgs = list("n_orderings" = 100))

setMethod("run", "BayesSpikeSolver",

  function (obj, target.gene, tfs, tf.weights=rep(1,length(tfs)), extraArgs=list()){

      # Check if target.gene is in the bottom 10% in mean expression; if so, send a warning      
      if(rowMeans(getAssayData(obj))[target.gene] < stats::quantile(rowMeans(getAssayData(obj)), probs = 0.1)){                   
          warning("Target gene mean expression is in the bottom 10% of all genes in the assay matrix")         
      }
            
      n_orderings <- 10
      
      # Check for n_orderings parameter and adopt a supplied one
      if("n_orderings" %in% names(extraArgs))
          n_orderings <- extraArgs[["n_orderings"]]
      
      # we don't try to handle tf self-regulation      
      deleters <- grep(target.gene, tfs)      
      if(length(deleters) > 0){          
          tfs <- tfs[-deleters]          
          tf.weights <- tf.weights[-deleters]          
          message(sprintf("BayesSpikeSolver removing target.gene from candidate regulators: %s", target.gene))          
      }      

    mtx <- getAssayData(obj)
    stopifnot(target.gene %in% rownames(mtx))
    stopifnot(all(tfs %in% rownames(mtx)))
    features <- t(mtx[tfs, ])
    target <- as.numeric(mtx[target.gene,])
    result <- vbsr(target, features, family='normal', n_orderings = n_orderings)
    tbl.out <- data.frame(beta=result$beta, pval=result$pval, z=result$z, post=result$post)
    rownames(tbl.out) <- tfs
    tbl.out$score <- -log10(tbl.out$pval)
    tbl.out <- tbl.out[order(tbl.out$score, decreasing=TRUE),]

    gene.cor <- sapply(rownames(tbl.out), function(tf) stats::cor(mtx[tf,], mtx[target.gene,]))
    tbl.out$gene.cor <- as.numeric(gene.cor)
    tbl.out
    })
#----------------------------------------------------------------------------------------------------