BASiCS_MCMC {BASiCS} | R Documentation

## BASiCS MCMC sampler

### Description

MCMC sampler to perform Bayesian inference for single-cell mRNA sequencing datasets using the model described in Vallejos et al (2015).

### Usage

    BASiCS_MCMC(
      Data,
      N,
      Thin,
      Burn,
      Regression,
      WithSpikes = TRUE,
      PriorParam = BASiCS_PriorParam(Data, PriorMu = "EmpiricalBayes"),
      SubsetBy = c("none", "gene", "cell"),
      NSubsets = 1,
      CombineMethod = c("pie", "consensus"),
      Weighting = c("naive", "n_weight", "inverse_variance"),
      Threads = getOption("Ncpus", default = 1L),
      BPPARAM = BiocParallel::bpparam(),
      ...
    )

### Arguments

`Data` | A `SingleCellExperiment` object. If `WithSpikes = TRUE`, this MUST be formatted to include the spike-ins and/or batch information (see vignette).
`N` | Total number of iterations for the MCMC sampler. Use `N>=max(4,Thin)`, `N` being a multiple of `Thin`.
`Thin` | Thining period for the MCMC sampler. Use `Thin>=2`.
`Burn` | Burn-in period for the MCMC sampler. Use `Burn>=1`, `Burn<N`, `Burn` being a multiple of `Thin`.
`Regression` | If `Regression = TRUE`, BASiCS exploits a joint prior formulation for mean and over-dispersion parameters to estimate a measure of residual over-dispersion is not confounded by mean expression. Recommended setting is `Regression = TRUE`.
`WithSpikes` | If `WithSpikes = TRUE`, BASiCS will use reads from added spike-ins to estimate technical variability. If `WithSpikess = FALSE`, BASiCS depends on replicated experiments (batches) to estimate technical variability. In this case, please supply the BatchInfo vector in `colData(Data)`. Default: `WithSpikes = TRUE`.
`PriorParam` | List of prior parameters for BASiCS_MCMC. Should be created using `BASiCS_PriorParam`.
`SubsetBy` | Character value specifying whether a divide and conquer inference strategy should be used. When this is set to `"gene"`, inference is performed on batches of genes separately, and when it is set to `"cell"`, inference is performed on batches of cells separately. Posterior distributions are combined using posterior interval estimation (see Li et al., 2016).
`NSubsets` | If `SubsetBy="gene"` or `SubsetBy="cell"`, `NSubsets` specifies the number of batches to create and perform divide and conquer inference with.
`CombineMethod` | The method used to combine subposteriors if `SubsetBy` is set to `"gene"` or `"cell"`. Options are `"pie"` corresponding to posterior interval estimation (see Li et al., 2016) or `"consensus"` (see Scott et al., 2016). Both of these methods use a form of weighted average to combine subposterior draws into the final posterior.
`Weighting` | The weighting method used in the weighted average chosen using `CombineMethod`. Available options are `"naive"` (unweighted), `"n_weight"` (weights are chosen based on the size of each partition) and `"inverse_variance"` (subposteriors are weighted based on the inverse of the variance of the subposterior for each parameter).
`Threads` | Integer specifying the number of threads to be used to parallelise parameter updates. Default value is the globally set `"Ncpus"` option, or 1 if this option is not set.
`BPPARAM` | A `BiocParallelParam` instance, used for divide and conquer inference.
`...` | Optional parameters. `AR` Optimal acceptance rate for adaptive Metropolis Hastings updates. It must be a positive number between 0 and 1. Default (and recommended): `AR = 0.44`. `StopAdapt` Iteration at which adaptive proposals are not longer adapted. Use `StopAdapt>=1`. Default: `StopAdapt = Burn`. `StoreChains` If `StoreChains = TRUE`, the generated `BASiCS_Chain` object is stored as a '.Rds' file (`RunName` argument used to index the file name). Default: `StoreChains = FALSE`. `StoreAdapt` If `StoreAdapt = TRUE`, trajectory of adaptive proposal variances (in log-scale) for all parameters is stored as a list in a '.Rds' file (`RunName` argument used to index file name). Default: `StoreAdapt = FALSE`. `StoreDir` Directory where output files are stored. Only required if `StoreChains = TRUE` and/or `StoreAdapt = TRUE`. Default: `StoreDir = getwd()`. `RunName` String used to index '.Rds' files storing chains and/or adaptive proposal variances. `PrintProgress` If `PrintProgress = FALSE`, console-based progress report is suppressed. `Start` Starting values for the MCMC sampler. We do not advise to use this argument. Default options have been tuned to facilitate convergence. If changed, it must be a list containing the following elements: `mu0`, `delta0`, `phi0`, `s0`, `nu0`, `theta0`, `ls.mu0`, `ls.delta0`, `ls.phi0`, `ls.nu0` and `ls.theta0` `GeneExponent/CellExponent` Exponents applied to the prior for MCMC updates. Intended for use only when performing divide & conquer MCMC strategies.

### Value

An object of class `BASiCS_Chain`.

### Author(s)

Catalina A. Vallejos cnvallej@uc.cl

Nils Eling eling@ebi.ac.uk

### References

Vallejos, Marioni and Richardson (2015). PLoS Computational Biology.

Vallejos, Richardson and Marioni (2016). Genome Biology.

Eling et al (2018). Cell Systems

Simple, Scalable and Accurate Posterior Interval Estimation Cheng Li and Sanvesh Srivastava and David B. Dunson arXiv (2016)

Bayes and Big Data: The Consensus Monte Carlo Algorithm Steven L. Scott, Alexander W. Blocker, Fernando V. Bonassi, Hugh A. Chipman, Edward I. George and Robert E. McCulloch International Journal of Management Science and Engineering Management (2016)

### Examples

    # Built-in simulated dataset
    set.seed(1)
    Data <- makeExampleBASiCS_Data()
    # To analyse real data, please refer to the instructions in:
    # https://github.com/catavallejos/BASiCS/wiki/2.-Input-preparation

    # Only a short run of the MCMC algorithm for illustration purposes
    # Longer runs migth be required to reach convergence
    Chain <- BASiCS_MCMC(Data, N = 50, Thin = 2, Burn = 10, Regression = FALSE,
                         PrintProgress = FALSE, WithSpikes = TRUE)

    # To run the regression version of BASiCS, use:
    Chain <- BASiCS_MCMC(Data, N = 50, Thin = 2, Burn = 10, Regression = TRUE,
                         PrintProgress = FALSE, WithSpikes = TRUE)

    # To run the non-spike version BASiCS requires the data to contain at least
    # 2 batches:
    set.seed(2)
    Data <- makeExampleBASiCS_Data(WithBatch = TRUE)
    Chain <- BASiCS_MCMC(Data, N = 50, Thin = 2, Burn = 10, Regression = TRUE,
                         PrintProgress = FALSE, WithSpikes = FALSE)

    # For illustration purposes we load a built-in 'BASiCS_Chain' object
    # (obtained using the 'BASiCS_MCMC' function)
    data(ChainSC)

    # `displayChainBASiCS` can be used to extract information from this output.
    # For example:
    head(displayChainBASiCS(ChainSC, Param = 'mu'))

    # Traceplot (examples only)
    plot(ChainSC, Param = 'mu', Gene = 1)
    plot(ChainSC, Param = 'phi', Cell = 1)
    plot(ChainSC, Param = 'theta', Batch = 1)

    # Calculating posterior medians and 95% HPD intervals
    ChainSummary <- Summary(ChainSC)

    # `displaySummaryBASiCS` can be used to extract information from this output
    # For example:
    head(displaySummaryBASiCS(ChainSummary, Param = 'mu'))

    # Graphical display of posterior medians and 95% HPD intervals
    # For example:
    plot(ChainSummary, Param = 'mu', main = 'All genes')
    plot(ChainSummary, Param = 'mu', Genes = 1:10, main = 'First 10 genes')
    plot(ChainSummary, Param = 'phi', main = 'All cells')
    plot(ChainSummary, Param = 'phi', Cells = 1:5, main = 'First 5 cells')
    plot(ChainSummary, Param = 'theta')

    # To constrast posterior medians of cell-specific parameters
    # For example:
    par(mfrow = c(1,2))
    plot(ChainSummary, Param = 'phi', Param2 = 's', SmoothPlot = FALSE)
    # Recommended for large numbers of cells
    plot(ChainSummary, Param = 'phi', Param2 = 's', SmoothPlot = TRUE)

    # To constrast posterior medians of gene-specific parameters
    par(mfrow = c(1,2))
    plot(ChainSummary, Param = 'mu', Param2 = 'delta', log = 'x',
         SmoothPlot = FALSE)
    # Recommended
    plot(ChainSummary, Param = 'mu', Param2 = 'delta', log = 'x',
         SmoothPlot = TRUE)

    # To obtain denoised rates / counts, see:
    # help(BASiCS_DenoisedRates)
    # and
    # help(BASiCS_DenoisedCounts)

    # For examples of differential analyses between 2 populations of cells see:
    # help(BASiCS_TestDE)

[Package BASiCS version 2.2.4 Index]