R/BASiCS_LoadChain.R

In catavallejos/BASiCS: Bayesian Analysis of Single-Cell Sequencing data

#' @title Loads pre-computed MCMC chains generated by the
#' \code{\link[BASiCS]{BASiCS_MCMC}} function
#'
#' @description  Loads pre-computed MCMC chains generated by the
#' \code{\link[BASiCS]{BASiCS_MCMC}} function, creating
#' a \code{\linkS4class{BASiCS_Chain}}  object
#'
#' @param RunName String used to index `.Rds` file containing the MCMC chain
#' (produced by the \code{\link[BASiCS]{BASiCS_MCMC}} function, with
#' \code{StoreChains = TRUE})
#' @param StoreDir Directory where `.Rds` file is stored.
#' Default: \code{StoreDir = getwd()}
#' @param StoreUpdatedChain Only required when the input files contain an
#' outdated version of a \code{\linkS4class{BASiCS_Chain}} object.
#' If \code{StoreUpdatedChain = TRUE}, an updated object is saved
#' (this overwrites original input file, if it was an `.Rds` file).
#'
#' @return An object of class \code{\linkS4class{BASiCS_Chain}}.
#'
#' @examples
#'
#' Data <- makeExampleBASiCS_Data()
#' Chain <- BASiCS_MCMC(
#'   Data,
#'   N = 50,
#'   Thin = 5,
#'   Burn = 5,
#'   Regression = FALSE,
#'   StoreChains = TRUE,
#'   StoreDir = tempdir(),
#'   RunName = "Test"
#' )
#' ChainLoad <- BASiCS_LoadChain(RunName = "Test", StoreDir = tempdir())
#' @seealso \code{\link[BASiCS]{BASiCS_Chain}}
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#'
#' @export
BASiCS_LoadChain <- function(RunName = "",
                             StoreDir = getwd(),
                             StoreUpdatedChain = FALSE) {
  
  if (file.exists(file.path(StoreDir, paste0("chain_", RunName, ".Rds")))) {
    Chain <- readRDS(file.path(StoreDir, paste0("chain_", RunName, ".Rds")))
    
    if (methods::.hasSlot(Chain, "mu")) {
      if (!is.null(Chain@mu)) {
        message(
          "`BASiCS_Chain` class definition was outdated. \n",
          "Object updated to be compatible with BASiCS version ",
          utils::packageVersion("BASiCS"), ".\n",
          "Set 'StoreUpdatedChain' = TRUE to save updated object.\n",
          "(this overwrites original input file).\n"
        )
        Chain <- .updateObject(Chain)
        if (StoreUpdatedChain) {
          saveRDS(
            Chain,
            file = file.path(StoreDir, paste0("chain_", RunName, ".Rds"))
          )
        }
      }
    }
    
  } else {
    file <- file.path(StoreDir, paste0("chain_mu_", RunName, ".txt"))
    if (file.exists(file)) {
      Mu <- read.delim(
        file,
        sep = " ",
        check.names = FALSE
      )
      rownames(Mu) <- NULL
      Delta <- read.delim(
        file.path(StoreDir, paste0("chain_delta_", RunName, ".txt")),
        sep = " ",
        check.names = FALSE
      )
      rownames(Delta) <- NULL
      Phi <- read.delim(
        file.path(StoreDir, paste0("chain_phi_", RunName, ".txt")),
        sep = " ",
        check.names = FALSE
      )
      rownames(Phi) <- NULL
      
      # Add-hoc fix for the no-spikes case
      file <- file.path(StoreDir, paste0("chain_s_", RunName, ".txt"))
      if (file.exists(file)) {
        S <- read.delim(
          file,
          sep = " ",
          check.names = FALSE
        )
        rownames(S) <- NULL        
      } else {
        S <- matrix(1, ncol = ncol(Phi), nrow = nrow(Phi))
      }
      
      Nu <- read.delim(
        file.path(StoreDir, paste0("chain_nu_", RunName, ".txt")),
        sep = " ",
        check.names = FALSE
      )
      rownames(Nu) <- NULL
      Theta <- read.delim(
        file.path(StoreDir, paste0("chain_theta_", RunName, ".txt")),
        sep = " ",
        check.names = FALSE
      )
      rownames(Theta) <- NULL
      
      Chain <- newBASiCS_Chain(
        list(
          "mu" = as.matrix(Mu),
          "delta" = as.matrix(Delta),
          "phi" = as.matrix(Phi),
          "s" = as.matrix(S),
          "nu" = as.matrix(Nu),
          "theta" = as.matrix(Theta)
        )
      )
      
      message(
        "`BASiCS_Chain` class definition was outdated. \n",
        "Object updated to be compatible with BASiCS version ",
        utils::packageVersion("BASiCS"), ".\n",
        "Set 'StoreUpdatedChain' = TRUE to save updated object.\n"
      )
      
      if (StoreUpdatedChain) {
        saveRDS(Chain,
                file = file.path(StoreDir, paste0("chain_", RunName, ".Rds"))
        )
      }
      
    } else {
      stop("Input file does not exist")
    }
  }
  return(Chain)
}