R/ARSyNbac.R

In ConesaLab/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. Use TRUE when the source of the batch effect is known.
#' @param filterNoise Logical. If TRUE (default) structured noise is removed form residuals. Use this option when there is an unknown source of batch effect in data.
#' @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 explained_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, filterNoise = TRUE,
#'                   Interaction=TRUE, Variability = 0.90, beta = 2,
#'                   modelName = "Model 1",
#'                   showplot = FALSE)
#'
ARSyNbac <- function(mbac, batchEstimation = TRUE, filterNoise = 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
    }
  } 
  if(filterNoise) {
    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)
  
}