R/Methods.R

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

# Remove functions that are not necessary anymore

#########################################################################
# Methods for BASiCS_Chain objects

#' @name BASiCS_Chain-methods
#' @aliases show,BASiCS_Chain-method
#'
#' @title 'show' method for BASiCS_Chain objects
#'
#' @description 'show' method for \code{\linkS4class{BASiCS_Chain}} objects.
#'
#' @param object A \code{\linkS4class{BASiCS_Chain}} object.
#'
#' @return Prints a summary of the properties of \code{object}.
#'
#' @examples
#'
#' Data <- makeExampleBASiCS_Data()
#' Chain <- BASiCS_MCMC(Data, N = 50, Thin = 2, Burn = 2, Regression = FALSE)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#'
#' @rdname BASiCS_Chain-methods
#' @export
setMethod(
  "show",
  signature = "BASiCS_Chain",
  definition = function(object) {

    N <- nrow(object@parameters$mu)
    q.bio <- ncol(object@parameters$delta)
    n <- ncol(object@parameters$nu)
    nBatch <- ncol(object@parameters$theta)
    if (!("parameters" %in% methods::slotNames(object))) {
      cat("'BASiCS_Chain' object was generated by an old version of BASiCS.\n")
    } else {
      N <- nrow(object@parameters$mu)
      q.bio <- ncol(object@parameters$delta)
      n <- ncol(object@parameters$nu)
      nBatch <- ncol(object@parameters$theta)
      cat("An object of class ", class(object), "\n",
          " ", N, " MCMC samples.\n", sep = "")
      if (nBatch > 1) {
        cat(" Dataset contains ", q.bio, " biological genes and ",
            n, " cells (", nBatch, " batches). \n", sep = "")
      } else {
        cat(" Dataset contains ", q.bio, " biological genes and ",
            n, " cells (1 batch). \n", sep = "")
      }
      ChainVersion <- object@.__classVersion__$BASiCS_Chain
      ChainVersion <- paste(ChainVersion, collapse = ".")
      cat(" Object stored using BASiCS version: ", ChainVersion, "\n",
          "Parameters: ", names(object@parameters), "\n")
    }
  }
)


.updateObject <- function(object, ..., verbose = FALSE) {
  if (is.null(object@mu)) {
    stop("Object was not created by an older version of BASiCS. \n")
  }

  New_Chain <- newBASiCS_Chain(
    parameters = list(
      mu = object@mu,
      delta = object@delta,
      phi = object@phi,
      s = object@s,
      nu = object@nu,
      theta = object@theta
    )
  )

  return(New_Chain)
}

#' @name BASiCS_Chain-methods
#' @aliases updateObject,BASiCS_Chain-method
#'
#' @title 'updateObject' method for BASiCS_Chain objects
#'
#' @description 'updateObject' method for \code{\linkS4class{BASiCS_Chain}}
#' objects. It is used to convert outdated \code{\linkS4class{BASiCS_Chain}}
#' objects into a version that is compatible with the Bioconductor release
#' of BASiCS. Do not use this method is \code{\linkS4class{BASiCS_Chain}}
#' already contains a \code{parameters} slot.
#'
#' @param ... Additional arguments of \code{\link[BiocGenerics]{updateObject}}
#' generic method. Not used within BASiCS.
#' @param verbose Additional argument of 
#' \code{\link[BiocGenerics]{updateObject}} generic method. Not used within 
#' BASiCS.
#'
#' @return Returns an updated \code{\linkS4class{BASiCS_Chain}} object that
#' contains all model parameters in a single slot object (list).
#'
#' @examples
#'
#' # Not run
#' # New_Chain <- updateObject(Old_Chain)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#'
#' @rdname BASiCS_Chain-methods
#' @export
setMethod(
  "updateObject",
  signature = "BASiCS_Chain",
  definition = .updateObject
)

#' @name Summary
#' @aliases Summary Summary,BASiCS_Chain-method
#'
#' @docType methods
#' @rdname Summary-BASiCS_Chain-method
#'
#' @title 'Summary' method for BASiCS_Chain objects
#'
#' @description For each of the BASiCS parameters (see Vallejos et al 2015),
#' \code{Summary} returns the corresponding postior medians and limits of
#' the high posterior density interval (probabilty equal to \code{prob})
#'
#' @param x A \code{\linkS4class{BASiCS_Chain}} object.
#' @param ... Unused, only included for consistency with the generic.
#' @param prob \code{prob} argument for \code{\link[coda]{HPDinterval}}
#' function.
#' @param na.rm Unused, only included for consistency with the generic.
#'
#' @return An object of class \code{\linkS4class{BASiCS_Summary}}.
#'
#' @examples
#'
#' data(ChainSC)
#' SummarySC <- Summary(ChainSC)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#' @importFrom stats na.omit
#' @export
setMethod("Summary",
          signature = "BASiCS_Chain",
          definition = function(x, ..., prob = 0.95, na.rm = FALSE) {
  param_ind <- !names(x@parameters) %in% c(
    "designMatrix",
    "RefFreq",
    "RBFLocations"
  )
  params <- x@parameters[param_ind]
  out <- lapply(params,
    function(n) {
      HPD <- matrix(
        data = NA,
        ncol = 3,
        nrow = ncol(n),
        dimnames = list(colnames(n), c("median", "lower", "upper"))
      )
      HPD[, 1] <- colMedians(n, na.rm = na.rm)
      ind_not_na <- !is.na(HPD[, 1])
      HPD[ind_not_na, 2:3] <- apply(n[, ind_not_na, drop=FALSE], 2, function(col) {
        if (na.rm) {
          col <- na.omit(col)
          # avoid coda error for no samples
          if (length(col) < 2) {
            return(NA)
          }
        }
        coda::HPDinterval(coda::mcmc(col), prob = prob)
      })
      HPD
    }
  )
  new("BASiCS_Summary", parameters = out)
})

#' @name subset
#' @aliases subset subset,BASiCS_Chain-method
#'
#' @docType methods
#' @rdname subset-BASiCS_Chain-method
#'
#' @title A 'subset' method for `BASiCS_Chain`` objects
#'
#' @description This can be used to extract a subset of a `BASiCS_Chain` object.
#' The subset can contain specific genes, cells or MCMC iterations
#'
#' @param x A \code{\linkS4class{BASiCS_Chain}} object.
#' @param Genes,Cells A vector of characters, logical values, or numbers, indicating
#'  which cells or genes will be extracted.
#' @param Iterations Numeric vector of positive integers indicating which MCMC
#' iterations will be extracted. The maximum value in \code{Iterations} must be
#' less or equal than the total number of iterations contained in the original
#' \code{\linkS4class{BASiCS_Chain}} object.
#'
#' @return An object of class \code{\linkS4class{BASiCS_Chain}}.
#'
#' @examples
#'
#' data(ChainSC)
#' 
#' # Extracts 3 first genes
#' ChainSC1 <- subset(ChainSC, Genes = rownames(ChainSC)[1:3])
#' # Extracts 3 first cells
#' ChainSC2 <- subset(ChainSC, Cells = colnames(ChainSC)[1:3])
#' # Extracts 10 first iterations
#' ChainSC3 <- subset(ChainSC, Iterations = 1:10)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#'
#' @export
setMethod(
  "subset",
  signature = "BASiCS_Chain",
  definition = function(x,
                        Genes = NULL,
                        Cells = NULL,
                        Iterations = NULL) {

    N <- nrow(displayChainBASiCS(x, Parameter = "mu"))
    q <- ncol(displayChainBASiCS(x, Parameter = "mu"))
    n <- ncol(displayChainBASiCS(x, Parameter = "nu"))
    GeneNames <- colnames(displayChainBASiCS(x, Parameter = "mu"))
    CellNames <- colnames(displayChainBASiCS(x, Parameter = "nu"))

    if (is.factor(Genes) | is.factor(Cells)) {
      stop("Factor index is not supported.")
    }
    if (is.numeric(Genes) | is.logical(Genes)) {
      if (is.logical(Genes) && length(Genes) != nrow(x)) {
        stop("Gene index length should match the number of features in the BASiCS_Chain object.")
      }
      Genes <- rownames(x)[Genes]
    }
    if (is.numeric(Cells) | is.logical(Cells)) {
      if (is.logical(Cells) && length(Cells) != nrow(x)) {
        stop("Cell index length should match the number of cells in the BASiCS_Chain object.")
      }
      Cells <- rownames(x)[Cells]
    }
    # Checking for valid arguments + assigning default values
    if (is.null(Iterations)) {
      Iterations <- seq_len(N)
    } else {
      stopifnot(all(Iterations == floor(Iterations)))
      stopifnot(min(Iterations) >= 1, max(Iterations) <= N )
    }
    if (is.null(Genes)) {
      Genes <- GeneNames
    } else {
      if (sum(!(Genes %in% GeneNames)) > 0) {
        stop("Some elements of 'Genes' are not present in the data")
      }
    }
    if (is.null(Cells)) {
      Cells <- CellNames
    } else {
      if (sum(!(Cells %in% CellNames)) > 0) {
        stop("Some elements of 'Cells' are not present in the data")
      }
    }

    # SelGenes <- GeneName %in% Genes
    # SelCells <- CellName %in% Cells
    Params <- names(x@parameters)

    out <- list()
    for(p in seq_len(length(Params))) {
      if (Params[p] %in% .OtherParams()) {
        out[[p]] <- get(Params[p], x@parameters)
      } else {
        if (Params[p] %in% c("theta", "beta", "sigma2")) {
          out[[p]] <- get(Params[p], x@parameters)[Iterations, , drop = FALSE]
        } else {
          if (Params[p] %in% c("mu", "delta", "epsilon")) {
            Sel <- Genes
          }
          if (Params[p] %in% c("phi", "s", "nu")) {
            Sel <- Cells
          }
          out[[p]] <- get(Params[p], x@parameters)[Iterations, Sel, drop = FALSE]
        }
      }
    }
    names(out) <- Params

    # Transform into BASiCS_Chain class object
    out1 <- newBASiCS_Chain(parameters = out)
    out1@.__classVersion__ <- x@.__classVersion__

    return(out1)
  }
)



#' @name dimnames
#' @aliases dimnames dimnames,BASiCS_Chain-method
#'
#' @docType methods
#' @rdname dimnames-BASiCS_Chain-method
#'
#' @title 'dimnames' method for BASiCS_Chain objects
#'
#' @description Returns the dimension names (genes x cells) of a 
#'  \code{\linkS4class{BASiCS_Chain}}
#'
#' @param x A \code{\linkS4class{BASiCS_Chain}} object.
#'
#' @return A list of two elements: (1) a vector of gene names and
#' (2) a vector of cell names.
#'
#' @examples
#'
#' data(ChainSC)
#' dimnames(ChainSC)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' 
#' @export
setMethod(
  "dimnames",
  signature = "BASiCS_Chain",
  definition = function(x) {
    list(colnames(x@parameters$mu), colnames(x@parameters$s))
  }
)


#' @name dim
#' @aliases dim dim,BASiCS_Chain-method
#'
#' @docType methods
#' @rdname dim-BASiCS_Chain-method
#'
#' @title 'dim' method for BASiCS_Chain objects
#'
#' @description Returns the dimensions (genes x cells) of a 
#'  \code{\linkS4class{BASiCS_Chain}}
#'
#' @param x A \code{\linkS4class{BASiCS_Chain}} object.
#'
#' @return An vector of dimensions
#'
#' @examples
#'
#' data(ChainSC)
#' dimnames(ChainSC)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' 
#' @export
setMethod(
  "dim",
  signature = "BASiCS_Chain",
  definition = function(x) c(ncol(x@parameters$mu), ncol(x@parameters$s))
)


# @importFrom ggExtra ggMarginal 
# @importFrom cowplot plot_grid
#' @name plot-BASiCS_Chain-method
#' @aliases plot,BASiCS_Chain-method plot,BASiCS_Chain,ANY-method
#'
#' @docType methods
#' @rdname plot-BASiCS_Chain-method
#'
#' @title 'plot' method for BASiCS_Chain objects
#'
#' @description 'plot' method for \code{\linkS4class{BASiCS_Chain}} objects
#'
#' @param x A \code{\linkS4class{BASiCS_Chain}} object.
#' @param Parameter Name of the slot to be used for the plot.
#' Possible values: \code{'mu'}, \code{'delta'},
#' \code{'phi'}, \code{'s'}, \code{'nu'}, \code{'theta'},
#' \code{'beta'}, \code{'sigma2'} and \code{'epsilon'}.
#' @param Gene Specifies which gene is requested.
#' Required only if \code{Parameter = 'mu'} or \code{'delta'}
#' @param Cell Specifies which cell is requested.
#' Required only if \code{Parameter = 'phi', 's'} or \code{'nu'}
#' @param Batch Specifies which batch is requested.
#' Required only if \code{Parameter = 'theta'}
#' @param RegressionTerm Specifies which regression coefficient is requested.
#' Required only if \code{Parameter = 'beta'}
#' @param ... Unused.
#'
#' @return A plot object
#'
#' @examples
#'
#' help(BASiCS_MCMC)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#' 
#' @export
setMethod(
  "plot",
  signature = "BASiCS_Chain",
  definition = function(x,
                        Parameter = "mu",
                        Gene = NULL,
                        Cell = NULL,
                        Batch = 1,
                        RegressionTerm = NULL,
                        ...) {

    if (Parameter %in% .GeneParams()) {
      if (is.null(Gene)) {
        stop("'Gene' value is required")
      }
      Column <- Gene
    } else if (Parameter %in% .CellParams()) {
      if (is.null(Cell)) {
        stop("'Cell' value is required")
      }
      Column <- Cell
    } else if (Parameter == "theta") {
      if (is.null(Batch)) {
        stop("'Batch' value is required")
      }
      Column <- Batch
    } else if (Parameter == "beta") {
      if (is.null(RegressionTerm)) {
        stop("'RegressionTerm' value is required")
      }
      Column <- RegressionTerm
    } else if (Parameter == "sigma2") {
      Column <- 1
    }
    DF1 <- data.frame(
      "Iteration" = seq_len(nrow(x@parameters[[Parameter]])),
      "Draws" = x@parameters[[Parameter]][, Column]
    )
    # Code inspired by https://tinyurl.com/vezqqud
    MyAcf <- acf(x@parameters[[Parameter]][,Column], plot = FALSE)
    DF2 <- with(MyAcf, data.frame(lag, acf))
    
    # Traceplot
    p1 <- ggplot(DF1) + 
      geom_point(
        aes(x = .data$Iteration, y = .data$Draws),
        col = grDevices::adjustcolor("white", 
        alpha.f = 0)
      ) + 
      geom_line(aes(x = .data$Iteration, y = .data$Draws)) + 
      labs(
        title = colnames(x@parameters[[Parameter]])[Column],
        x = "Iteration", 
        y = "Parameter value"
      ) + 
      theme_classic()
    p1 <- ggExtra::ggMarginal(p1, type = "histogram", margins = "y")
      
    p2 <- ggplot(DF2, aes(x = .data$lag, y = .data$acf)) +
      theme_classic() + 
      geom_hline(aes(yintercept = 0)) + 
      geom_segment(
        mapping = aes(xend = .data$lag, yend = 0)
      ) +
      ggtitle(colnames(x@parameters[[Parameter]])[Column])
    
    cowplot::plot_grid(p1, p2)
            
  }
)

#' @name displayChainBASiCS-BASiCS_Chain-method
#' @aliases displayChainBASiCS displayChainBASiCS,BASiCS_Chain-method
#'
#' @docType methods
#' @rdname displayChainBASiCS-BASiCS_Chain-method
#'
#' @title Accessors for the slots of a BASiCS_Chain object
#'
#' @description Accessors for the slots of a \code{\linkS4class{BASiCS_Chain}}
#'
#' @param object an object of class \code{\linkS4class{BASiCS_Chain}}
#' @param Parameter Name of the slot to be used for the accessed.
#' Possible values: \code{'mu'}, \code{'delta'}, \code{'phi'},
#' \code{'s'}, \code{'nu'}, \code{'theta'}, \code{'beta'},
#' \code{'sigma2'} and \code{'epsilon'}.
#'
#' @return The requested slot of a \code{\linkS4class{BASiCS_Chain}} object
#'
#' @examples
#'
#' help(BASiCS_MCMC)
#'
#' @seealso \code{\linkS4class{BASiCS_Chain}}
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#'
#' @export
setMethod(
  "displayChainBASiCS",
  signature = "BASiCS_Chain",
  definition = .GetParam
)

################################################################################
# Methods for BASiCS_Summary objects

#' @name BASiCS_Summary-methods
#' @aliases show,BASiCS_Summary-method
#'
#' @title 'show' method for BASiCS_Summary objects
#'
#' @description 'show' method for \code{\linkS4class{BASiCS_Summary}} objects.
#'
#' @param object A \code{\linkS4class{BASiCS_Summary}} object.
#'
#' @return Prints a summary of the properties of \code{object}.
#'
#' @examples
#'
#' data(ChainSC)
#' show(ChainSC)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#'
#' @rdname BASiCS_Summary-methods
#' 
#' @export
setMethod(
  "show",
  signature = "BASiCS_Summary",
  definition = function(object) {
    q <- nrow(object@parameters$mu)
    q.bio <- nrow(object@parameters$delta)
    n <- nrow(object@parameters$phi)
    nBatch <- nrow(object@parameters$theta)
    cat("An object of class", class(object), "\n",
        " Contains posterior medians and the limits of the \n",
        " HPD interval for all BASiCS model parameters.\n")

    if (nBatch > 1) {
      cat("  Dataset contains ", q, " genes ",
          "(", q.bio, " biological and ", q-q.bio, " technical) \n",
          "  and ", n, " cells (", nBatch, " batches). \n", sep = "")
    }
    else {
      cat("  Dataset contains ", q, " genes ",
          "(", q.bio, " biological and ", q-q.bio, " technical) \n",
          "  and ", n, " cells (1 batch). \n", sep = "")
    }
    ChainVersion <- object@.__classVersion__$BASiCS_Summary
    ChainVersion <- paste0(ChainVersion[1], ".",
                           ChainVersion[2], ".",
                           ChainVersion[3])
    cat("  Object stored using BASiCS version: ", ChainVersion, "\n",
        " Parameters: ", names(object@parameters), "\n")
  }
)




#' @name plot-BASiCS_Summary-method
#' @aliases plot,BASiCS_Summary-method, plot,BASiCS_Summary,ANY-method
#'
#' @docType methods
#' @rdname plot-BASiCS_Summary-method
#'
#' @title 'plot' method for BASiCS_Summary objects
#'
#' @description 'plot' method for \code{\linkS4class{BASiCS_Summary}} objects
#'
#' @param x A \code{\linkS4class{BASiCS_Summary}} object.
#' @param Param Name of the slot to be used for the plot.
#' Possible values: \code{'mu'}, \code{'delta'},
#' \code{'phi'}, \code{'s'}, \code{'nu'}, \code{'theta'}, \code{'beta'},
#' \code{'sigma2'} and \code{'epsilon'}.
#' @param Param2 Name of the second slot to be used for the plot.
#' Possible values: \code{'mu'}, \code{'delta'}, \code{'epsilon'},
#' \code{'phi'}, \code{'s'} and \code{'nu'}
#' (combinations between gene-specific and
#' cell-specific parameters are not admitted).
#' @param Genes Specifies which genes are requested.
#' Required only if \code{Param = 'mu'}, \code{'delta'} or \code{'epsilon'}.
#' @param Cells Specifies which cells are requested.
#' Required only if \code{Param = 'phi', 's'} or \code{'nu'}
#' @param Batches Specifies which batches are requested.
#' Required only if \code{Param = 'theta'}
#' @param RegressionTerms Specifies which regression coefficients are requested.
#' Required only if \code{Param = 'beta'}
#'
#' @param xlab As in \code{\link[graphics]{par}}.
#' @param ylab As in \code{\link[graphics]{par}}.
#' @param xlim As in \code{\link[graphics]{par}}.
#' @param ylim As in \code{\link[graphics]{par}}.
#' @param pch As in \code{\link[graphics]{par}}.
#' @param col As in \code{\link[graphics]{par}}.
#' @param bty As in \code{\link[graphics]{par}}.
#' @param SmoothPlot Logical parameter. If \code{TRUE},
#' transparency will be added to the color of the dots.
#' @param ... Other graphical parameters (see \code{\link[graphics]{par}}).
#'
#' @return A plot object
#'
#' @examples
#'
#' help(BASiCS_MCMC)
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#'
#' @export
setMethod(
  "plot",
  signature = "BASiCS_Summary",
  definition = function(x,
                        Param = "mu",
                        Param2 = NULL,
                        Genes = NULL,
                        Cells = NULL,
                        Batches = NULL,
                        RegressionTerms = NULL,
                        xlab = "",
                        ylab = "",
                        xlim = "",
                        ylim = NULL,
                        pch = 16,
                        col = "blue",
                        bty = "n",
                        SmoothPlot = TRUE,
                        ...) {

    if (!(Param %in% names(x@parameters))) {
      stop("'Param' argument is invalid")
    }

    q <- nrow(x@parameters$mu)
    q.bio <- nrow(x@parameters$delta)
    n <- nrow(x@parameters$s)
    nBatch <- nrow(x@parameters$theta)

    if (is.null(Genes)) {
      Genes <- seq_len(q.bio)
    }
    if (is.null(Cells)) {
      Cells <- seq_len(n)
    }
    if (is.null(Batches)) {
      Batches <- seq_len(nBatch)
    }

    if("beta" %in% names(x@parameters)) {
      k <- nrow(x@parameters$beta)
      if (is.null(RegressionTerms)) {
        RegressionTerms <- seq_len(k)
      }
    }

    if (is.null(Param2)) {
      object <- x@parameters[[Param]]

      if (Param %in% .GeneParams()) {
        Columns <- Genes
        ylabInd <- "i"
        if (xlab == "") {
          xlab <- "Gene"
        }
      } else if (Param %in% .CellParams()) {
        Columns <- Cells
        ylabInd <- "j"
        if (xlab == "") {
          xlab <- "Cell"
        }
      } else if (Param == "theta") {
        Columns <- Batches
        ylabInd <- "b"
        if (xlab == "") {
          xlab <- "Batch"
        }
      } else if (Param == "beta") {
        Columns <- RegressionTerms
        ylabInd <- "k"
        if (xlab == "") {
          xlab <- "Regression term"
        }
      } else if (Param == "sigma2") {
        object <- x@parameters$sigma2
        Columns <- 1
        ylabInd <- "i"
        if (xlab == "") {
          xlab <- "Gene"
        }
      }

      if (ylab == "") {
        ylab <- bquote(.(Param)[.(ylabInd)])
      }
      if (is.null(ylim)) {
        ylim = c(min(object[Columns, 2]), max(object[Columns, 3]))
      }

      # Point estimates
      plot(x = Columns,
           y = object[Columns, 1],
           xlab = xlab,
           ylab = ylab,
           ylim = ylim,
           pch = pch,
           col = col,
           bty = bty,
           ...
      )

      # Add HPD interval
      for (Column in Columns) {
        BarLength <- ifelse(length(Columns) <= 10,
                            0.1,
                            2 / length(Columns))

        segments(x0 = Column,
                 y0 = object[Column, 2],
                 y1 = object[Column, 3],
                 col = col)
        segments(x0 = Column - BarLength,
                 y0 = object[Column, 2],
                 x1 = Column + BarLength,
                 col = col)
        segments(x0 = Column - BarLength,
                 y0 = object[Column, 3],
                 x1 = Column + BarLength,
                 col = col)
      }

    } else {
      if (!Param2 %in% names(x@parameters)) {
        stop("'Param2' is not a parameter in this BASiCS_Summary object")
      }
      .CheckValidCombination(Param, Param2)
      if (SmoothPlot) {
        col <- grDevices::rgb(grDevices::col2rgb(col)[1],
                              grDevices::col2rgb(col)[2],
                              grDevices::col2rgb(col)[3],
                              50,
                              maxColorValue = 255)
      }
      if (Param %in% .GeneParams()) {
        Columns <- Genes
        ylabInd <- "i"
      }
      if (Param %in% .CellParams()) {
        Columns <- Cells
        ylabInd <- "j"
      }
      graphics::plot(
        x@parameters[[Param]][Columns, 1],
        x@parameters[[Param2]][Columns, 1],
        xlab = bquote(.(Param)[.(ylabInd)]),
        ylab = bquote(.(Param2)[.(ylabInd)]),
        xlim = c(
          min(x@parameters[[Param]][Columns, 1]),
          max(x@parameters[[Param]][Columns, 1])
        ),
        ylim = c(
          min(x@parameters[[Param2]][Columns, 1]),
          max(x@parameters[[Param2]][Columns, 1])
        ),
        pch = pch,
        col = col,
        bty = bty,
        ...
      )
    }
  }
)


#' @name displaySummaryBASiCS-BASiCS_Summary-method
#' @aliases displaySummaryBASiCS displaySummaryBASiCS,BASiCS_Summary-method
#'
#' @docType methods
#' @rdname displaySummaryBASiCS-BASiCS_Summary-method
#'
#' @title Accessors for the slots of a \code{\linkS4class{BASiCS_Summary}}
#' object
#'
#' @description Accessors for the slots of a \code{\linkS4class{BASiCS_Summary}}
#' object
#'
#' @param object an object of class \code{\linkS4class{BASiCS_Summary}}
#' @param Parameter Name of the slot to be used for the accessed.
#' Possible values: \code{'mu'}, \code{'delta'}, \code{'phi'},
#' \code{'s'}, \code{'nu'}, \code{'theta'}, \code{'beta'},
#' \code{'sigma2'} and \code{'epsilon'}.
#'
#' @return The requested slot of a \code{\linkS4class{BASiCS_Summary}} object
#'
#' @examples
#'
#' help(BASiCS_MCMC)
#'
#' @seealso \code{\linkS4class{BASiCS_Summary}}
#'
#' @author Catalina A. Vallejos \email{cnvallej@@uc.cl}
#' @author Nils Eling \email{eling@@ebi.ac.uk}
#'
#' @include utils_Misc.R
#' @export
setMethod(
  "displaySummaryBASiCS",
  signature = "BASiCS_Summary",
  definition = .GetParam
)


#' @docType methods
#' @rdname show-BASiCS_ResultsDE-method
#'
#' @title Accessors for the slots of a \code{\linkS4class{BASiCS_ResultsDE}}
#' object
#'
#' @param object an object of class \code{\linkS4class{BASiCS_ResultsDE}}
#' @return Prints a summary of the properties of \code{object}.
#' @examples
#'
#' help(BASiCS_MCMC)
#'
#' @seealso \code{\link[methods]{show}}
#' @export
setMethod("show", signature = "BASiCS_ResultsDE", 
  function(object) {
    cat("An object of class BASiCS_ResultsDE, containing:\n")
    for (x in object@Results) {
      print(x)
    }
  }
)

#' @docType methods
#' @rdname show-BASiCS_ResultVG-method
#'
#' @title Accessors for the slots of a \code{\linkS4class{BASiCS_ResultVG}}
#' object
#'
#' @param object an object of class \code{\linkS4class{BASiCS_ResultsDE}}
#' @return Prints a summary of the properties of \code{object}.
#' @examples
#'
#' help(BASiCS_MCMC)
#'
#' @seealso \code{\link[methods]{show}}
#' @export
setMethod("show", signature = "BASiCS_ResultVG", 
  function(object) {
    cat("An object of class BASiCS_ResultVG.\n")
    if (object@Method == "Percentile") {
      cat(
        sum(.VG(object), na.rm = TRUE), 
        "genes classified as", object@Name, "using: \n",
        "- The", round(100 * object@Threshold, 2), 
        "percentile of variable genes \n",
        "- Evidence threshold =", object@ProbThreshold, "\n",
        "- EFDR =", round(100 * object@EFDR, 2), "% \n",
        "- EFNR =", round(100 * object@EFNR, 2), "% \n"
      )
    } else if (object@Method == "Variance") {
      cat(
        sum(.VG(object), na.rm = TRUE),
        " genes classified as", object@Name, "using:\n",
        "- Variance contribution threshold =",
        round(100 * object@Threshold, 2), "% \n",
        "- Evidence threshold =", object@ProbThreshold, "\n",
        "- EFDR =", round(100 * object@EFDR, 2), "% \n",
        "- EFNR =", round(100 * object@EFNR, 2), "% \n"
      )
    }
  }
)


#' @docType methods
#' @rdname show-BASiCS_ResultDE-method
#'
#' @title Accessors for the slots of a \code{\linkS4class{BASiCS_ResultDE}}
#' object
#'
#' @param object an object of class \code{\linkS4class{BASiCS_ResultDE}}
#' @return Prints a summary of the properties of \code{object}.
#'
#' @examples
#'
#' help(BASiCS_MCMC)
#'
#' @seealso \code{\link[methods]{show}}
#' @export
setMethod("show", signature = "BASiCS_ResultDE", 
  function(object) {
    cat("-------------------------------------------------------------\n")
    cat("  An object of class BASiCS_ResultDE.\n")
    # measures <- c("Mean", "Disp", "ResDisp")
    diffName <- paste0("ResultDiff", object@Name)
    
    nPlus1 <- sum(object@Table[[diffName]] == paste0(object@GroupLabel1, "+"))
    nPlus2 <- sum(object@Table[[diffName]] == paste0(object@GroupLabel2, "+"))
    message <- paste(
      "-------------------------------------------------------------\n",
      nPlus1 + nPlus2, "genes with a change in",
      .MeasureName(object@Name), "\n",
      "- Higher", .MeasureName(object@Name), "in", 
      object@GroupLabel1, "samples:", nPlus1, "\n",
      "- Higher", .MeasureName(object@Name), "in", 
      object@GroupLabel2, "samples:", nPlus2, "\n",
      "-", .cap(.DistanceName(object@Name)), 
      "tolerance =", round(2^(object@Epsilon) * 100, 2), "% \n",
      "- Probability threshold =", object@ProbThreshold, "\n",
      "- EFDR =", round(100 * object@EFDR, 2), "% \n",
      "- EFNR =", round(100 * object@EFNR, 2), "%"
    )
    if (object@Name == "delta") {
      NotDE <- sum(object@Table[[diffName]]  == "ExcludedFromTesting")
      message <- paste(
        message, "\n",
        "NOTE: differential dispersion assessment only applied to the \n",
        sum(NotDE), "genes for which the mean did not change. \n",
        "and that were included for testing."
      )
    } else if (object@Name == "epsilon") {
      NotExcluded <- sum(object@Table[[diffName]] != "ExcludedFromTesting")
      sn <- sum(NotExcluded)
      message <- paste(
        message, "\n",
        "NOTE: differential residual dispersion assessment applied to\n",
        sn, "genes expressed in at least 2 cells per condition\n",
        "and that were included for testing."
      )
    }
    cat(message, "\n")
  }
)

#' @rdname subset-methods
#' @export
setMethod("[",
  signature = signature("BASiCS_ResultsDE", "ANY", "ANY", "ANY"),
  function(x, i, j, drop = FALSE) {
    if (missing(i)) {
      i <- seq_len(nrow(x@Results[[1]]@Table))
    } else if (!is.character(i)) {
      stop("Only GeneName values should be used as row indices")
    }
    if (missing(j)) {
      j <- seq_len(ncol(x@Results[[1]]@Table))
    }
    ir <- match(i, x@RowData$GeneName)
    x@RowData <- x@RowData[ir, seq_len(ncol(x@RowData)), drop = FALSE]
    x@Results <- lapply(
      x@Results,
      `[`,
      i = i,
      j = j,
      drop = drop
    )
    x
  }
)

#' Methods for subsetting \linkS4class{BASiCS_Result} and 
#' \linkS4class{BASiCS_ResultsDE} objects.
#' @param x Object being subsetted.
#' @param i See \code{?`[`}, \code{?`[[`}
#' @param j,drop Ignored.
#' @return An object of the same class as \code{x}.
#' @rdname subset-methods
#' @export
setMethod("[[",
  signature = signature("BASiCS_ResultsDE", "ANY"),
  function(x, i) {
    x@Results[[i]]
  }
)

#' @rdname subset-methods
#' @export
setMethod("[",
  signature = signature("BASiCS_Result", "ANY", "ANY", "ANY"),
  function(x, i, j, drop = FALSE) {
    if (missing(i)) {
      i <- seq_len(nrow(x@Table))
    } else if (is.character(i)) {
      i <- match(i, x@Table$GeneName)
    }
    j <- seq_len(ncol(x@Table))
    x@Table <- x@Table[i, j, drop = FALSE]
    x
  }
)


#' Methods for formatting \linkS4class{BASiCS_Result} and 
#' \linkS4class{BASiCS_ResultsDE} objects.
#' @param x Object being subsetted.
#' @param Parameter Character scalar indicating which of the 
#' \linkS4class{BASiCS_Result} should be formatted.
#' @param Filter Logical scalar indicating whether results should be
#' filtered based on differential expression or HVG/LVG status if
#' \code{ProbThreshold=NULL}, or a probability threshold if
#' \code{ProbThreshold=NULL}
#' @param ProbThreshold Probability threshold to be used if \code{Filter=TRUE}
#' @param ... Passed to \code{\link{format}}.
#' @return A \code{data.frame}.
#' @rdname format-methods
#' @export
setMethod("format",
  signature = signature("BASiCS_ResultsDE"),
  function(x, Parameter, Filter = TRUE, ProbThreshold = NULL, ...) {
    .Deprecated("as.data.frame", old = "format")
    format(
      as.data.frame(
        x,
        Parameter = Parameter,
        Filter = Filter,
        ProbThreshold = ProbThreshold
      ),
      ...
    )
  }
)



#' Converting BASiCS results objects to data.frames
#' 
#' @param x An object of class \linkS4class{BASiCS_ResultVG}, 
#'  \linkS4class{BASiCS_ResultDE}, or \linkS4class{BASiCS_ResultsDE}.
#' @param Parameter For \linkS4class{BASiCS_ResultsDE} objects only. 
#'  Character scalar specifying which table of results to 
#'  output. Available options are "Mean", (mu, mean expression),
#'  "Disp" (delta, overdispersion) and "ResDisp" 
#'  (epsilon, residual overdispersion).
#' @param Filter Logical scalar. If \code{TRUE}, output only entries
#' corresponding to genes that pass the decision rule used in the probabilistic
#' test.
#' @param ProbThreshold Only used if \code{Filter=TRUE}.
#' Numeric scalar specifying the probability threshold to be used when filtering
#' genes. Default is to use the threshold used in the original decision rule
#' when the test was performed.
#' @return A data.frame of test results.
#' @rdname as.data.frame-methods
#' @importFrom BiocGenerics as.data.frame
#' @export
setMethod("as.data.frame", signature = signature("BASiCS_ResultsDE"),
  function(x, Parameter, Filter = TRUE, ProbThreshold = NULL) {
    Parameter <- match.arg(Parameter, choices = names(x@Results))
    merge(
      as.data.frame(
        x@Results[[Parameter]],
        Filter = Filter,
        ProbThreshold = ProbThreshold
      ),
      x@RowData,
      by = "GeneName"
    )
  }
)


#' @rdname format-methods
#' @export
setMethod("format",
  signature = signature("BASiCS_ResultDE"),
  function(x, Filter = TRUE, ProbThreshold = NULL, ...) {
    .Deprecated("as.data.frame", old = "format")
    format(
      as.data.frame(x, Filter = Filter, ProbThreshold = ProbThreshold),
      ...
    )
  }
)


#' @rdname as.data.frame-methods
#' @export
setMethod("as.data.frame", signature = signature("BASiCS_ResultDE"),
  function(x, Filter = TRUE, ProbThreshold = NULL) {
    if (Filter) {
      if (is.null(ProbThreshold)) {
        ProbThreshold <- x@ProbThreshold
      }
      ind <- x@Table$Prob > ProbThreshold & .WhichDiffRes(x)
    } else {
      ind <- seq_len(nrow(x@Table))
    }
    x@Table[ind, ]
  }
)

#' @rdname format-methods
#' @export
setMethod("format",
  signature = signature("BASiCS_ResultVG"),
  function(x, Filter = TRUE, ProbThreshold = NULL, ...) {
    .Deprecated("as.data.frame", old = "format")
    format(
      as.data.frame(x, Filter = Filter, ProbThreshold = ProbThreshold),
      ...
    )
  }
)

#' @rdname as.data.frame-methods
#' @export
setMethod("as.data.frame", signature = signature("BASiCS_ResultVG"),
  function(x, Filter = TRUE, ProbThreshold = NULL) {
    if (Filter) {
      if (is.null(ProbThreshold)) {
        ind <- .VG(x)
      } else {
        ind <- x@Table$Prob > ProbThreshold
      }
    } else {
      ind <- seq_len(nrow(x@Table))
    }
    merge(
      x@Table[ind, ],
      x@RowData,
      by = "GeneName"
    )
  }
)

#' rowData getter and setter for \linkS4class{BASiCS_ResultsDE} 
#' and \linkS4class{BASiCS_ResultVG} objects.
#' @param x \linkS4class{BASiCS_ResultVG} or \linkS4class{BASiCS_ResultsDE}
#' object.
#' @param value New \code{rowData} value for setter method.
#' @return For the getter, a \linkS4class{DFrame}. For setter, the modified
#' \code{x}.
#' @rdname rowData
#' @export
setMethod("rowData",
  signature = signature("BASiCS_ResultsDE"),
  function(x) {
    x@RowData
  }
)

#' @rdname rowData
#' @export
setMethod("rowData<-",
  signature = signature("BASiCS_ResultsDE"),
  function(x, value) {
    x@RowData <- value
    x
  }
)

#' @rdname rowData
#' @export
setMethod("rowData",
  signature = signature("BASiCS_ResultVG"),
  function(x) {
    x@RowData
  }
)

#' @rdname rowData
#' @export
setMethod("rowData<-",
  signature = signature("BASiCS_ResultVG"),
  function(x, value) {
    x@RowData <- value
    x
  }
)