R/ARSyNbac.R

In MultiBaC: Multiomic Batch effect Correction

#' @include ASCA2f.R
#' @include ASCA1f.R
#' @include ARSyNcomponents.R
#' @include auxfunctions.R
#'
NULL

#' ARSyNbac
#'
#' @param mbac mbac object generated by *createMbac*.
#' @param Interaction Logical. Whether to model the interaction between factors or not (FALSE by default).
#' @param Variability From 0 to 1. Minimum percent of data variability that must be explained by each model. By default, 0.90.
#' @param batchEstimation Logical. If TRUE (default) the batch effect is estimated and used to correct the data. If batch effect is unknown or it is not the main source of noise, this argument must be set to FALSE and ARSyNbac will extract unwanted effects from residuals.
#' @param beta Numeric. Components that represent more than beta times the average variability are identified as systematic noise in residuals. Used in noise reduction mode. By default, 2.
#' @param modelName Name of the model created. This name will be showed if you use the explaine_var plot function. By default, "Model 1".
#' @param showplot Logical. If TRUE (default), the explained_var plot is showed. This plot represents the number of components selected for the ARSyN model.
#'
#' @return Custom mbac object. Elements in a mbac object:
#' \enumerate{
#'     \item ListOfBatches: A list of MultiAssayExperiment objects (one per batch).
#'     \item commonOmic Name of the common omic between the batches. It must be one of the names in omicNames argument. If NULL (default), the omic names that appears more times is selected as commonOmic.
#'     \item CorrectedData: Same structure than ListOfBatches but with the corrected data instead of the original.
#'     \item ARSyNmodels: ARSyN models created during MultiBaC performance (one per omic data type).
#' }
#'
#' @references
#' Nueda MJ, Ferrer A, Conesa A. ARSyN: A method for the identification and removal of systematic noise in multifactorial time course microarray experiments. Biostatistics. 2012;13:553–66.
#'
#' @export
#'
#' @examples
#' data('multiyeast')
#'
#' my_mbac <- createMbac (inputOmics = list(A.rna, B.rna, C.rna),
#'                        batchFactor = c("A", "B", "C"),
#'                        experimentalDesign = list("A" =  c("Glu+",
#'                        "Glu+", "Glu+", "Glu-",
#'                        "Glu-", "Glu-"),
#'                        "B" = c("Glu+", "Glu+", "Glu-", "Glu-"),
#'                        "C" = c("Glu+", "Glu+", "Glu-", "Glu-")),
#'                        omicNames = "RNA")
#' my_final_mbac <- ARSyNbac (my_mbac, batchEstimation = TRUE,
#'                   Interaction=TRUE, Variability = 0.90, beta = 2,
#'                   modelName = "Model 1",
#'                   showplot = FALSE)
#'
ARSyNbac <- function(mbac, batchEstimation = TRUE,
                     Interaction=FALSE, Variability = 0.90, beta = 2,
                     modelName = "Model 1",
                     showplot = TRUE) {

  ####################################
  ### --- Compute Inputs for ASCA
  ####################################

  # Input data
  X <- NULL
  experimentalDesign <- NULL
  batch_factor <- c()
  for (l in seq_along(mbac$ListOfBatches)) {
    X <- rbind(X, t(mbac$ListOfBatches[[l]]@ExperimentList[[1]]))
    experimentalDesign <- rbind(experimentalDesign, as.data.frame(mbac$ListOfBatches[[l]]@colData))
    batch_factor <- c(batch_factor, rep(l, dim(mbac$ListOfBatches[[l]][[1]])[2]))
  }

  # Experimental design factor
  designA <-  stats::model.matrix (~ 0 + factor(experimentalDesign[,1]))

  # Batch factor
  if (batchEstimation) {
    designB <- stats::model.matrix(~ 0 + factor(batch_factor))
  }

  if (batchEstimation) {
    Fac0<-c(1,2,2)
    if (Interaction) {
      names(Fac0) <- c("Model.a", "Model.bab", "Model.res")
    } else {
      names(Fac0) <- c("Model.a", "Model.b", "Model.res")
    }
    asca0<- ASCA.2f(X=X,
                    Designa=designA,
                    Designb=designB,
                    Fac=Fac0,
                    Interaction=Interaction,
                    Variability = Variability,
                    beta = beta)
    Fac<- ARSyNcomponents(asca0,Variability=Variability, beta = beta)
    for (i in 1:length(Fac)){
      Fac0[names(Fac[i])]<-Fac[names(Fac[i])]
    }
    asca<- ASCA.2f(X=X,
                   Designa=designA,
                   Designb=designB,
                   Fac=Fac0,Interaction=Interaction,
                   Variability = Variability,
                   beta = beta)
  } else {
    Fac0<-c(1,2)
    names(Fac0)<-c("Model.a","Model.res")
    asca0<- ASCA.1f(X=X,
                    Designa=designA,
                     Fac=Fac0,
                    Variability = Variability,
                    beta = beta)
    Fac<- ARSyNcomponents(asca0,Variability=Variability,beta=beta)
    for (i in 1:length(Fac)){
      Fac0[names(Fac[i])]<-Fac[names(Fac[i])]
    }
    asca<- ASCA.1f(X=X,
                   Designa=designA,
                   Fac=Fac0,
                   Variability = Variability,
                   beta = beta)
  }

  ####################################
  ### --- Writing filtered matrix
  ####################################

  X.filtered <- X

  if (batchEstimation) {
    if ( Interaction ) {
      X.filtered <- X.filtered - asca$Model.bab$TP # - asca$Model.b$TP
    } else {
      X.filtered <- X.filtered - asca$Model.b$TP
    }
  } else {
    X.filtered <- X.filtered - asca$Model.res$TP
  }

  data.filtered <- t(X.filtered)
  asca <- list(asca); names(asca) <- modelName

  # Separate matrix into original batches

  # output object
  output <- mbac$ListOfBatches
  roff <- 0
  for (lab in names(mbac$ListOfBatches)) {
    output[[lab]]@ExperimentList[[1]] <- data.filtered[,(roff+1):(roff + dim(output[[lab]]@ExperimentList[[1]])[2])]
    roff <- roff + dim(output[[lab]]@ExperimentList[[1]])[2]
  }

  retobj <- list("ListOfBatches" = mbac$ListOfBatches,
                 "CorrectedData" = c(lapply(output, function(x) x)),
                 "ARSyNmodels" = asca)

  mbac_out <- mbacClass(retobj)

  if(showplot) {
    plot.mbac(mbac_out, plotType = "var")
  }

  return(mbac_out)

}