BASiCS 1.2.0

Single-cell mRNA sequencing can uncover novel cell-to-cell heterogeneity in gene expression levels within seemingly homogeneous populations of cells. However, these experiments are prone to high levels of technical noise, creating new challenges for identifying genes that show genuine heterogeneous expression within the group of cells under study.

BASiCS (**B**ayesian **A**nalysis of **Si**ngle-**C**ell **S**equencing data)
is an integrated Bayesian hierarchical model that propagates statistical
uncertainty by simultaneously performing data normalisation (global scaling),
technical noise quantification and two types of **supervised** downstream
analyses:

**For a single group of cells**(Vallejos, Marioni, and Richardson 2015): BASiCS provides a criterion to identify highly (and lowly) variable genes within the group.**For two (or more) groups of cells**: BASiCS allows the identification of changes in gene expression between the groups. As in traditional differential expression tools, BASiCS can uncover changes in mean expression between the groups. Besides this, BASiCS can also uncover changes in gene expression variability in terms of:

**Over-dispersion**(Vallejos, Richardson, and Marioni 2016) — a measure for the excess of cell-to-cell variability that is observed with respect to Poisson sampling, after accounting for technical noise. This feature has led, for example, to novel insights in the context of immune cells across aging (Martinez-Jimenez et al. 2017). However, due to the strong mean/over-dispersion confounding that is typically observed for scRNA-seq datasets, the assessment of changes in over-dispersion is restricted to genes for which mean expression does not change between the groups.**Residual over-dispersion**(Eling et al. 2017) — a residual measure of variability given by departures with respect to a global mean/over-dispersion trend. Positive residual over-dispersion indicates that a gene exhibits more variation than expected relative to genes with similar expression levels; negative residual over-dispersion suggests less variation than expected. To use this feature, please set`Regression = TRUE`

as a function parameter in`BASiCS_MCMC`

.

In all cases, a probabilistic output is provided and a decision rule is calibrated using the expected false discovery rate (EFDR) (Newton et al. 2004).

A brief description for the statistical model implemented in BASiCS is
provided in Section 6 of this document. The original
implementation of BASiCS (Vallejos, Marioni, and Richardson 2015) requires the use of **spike-in**
molecules — that are artificially introduced to each cell’s lysate
— to perform these analyses. More recently, Eling et al. (2017) extendeded BASiCS
to also address datasets for which spikes-ins are not available
(see Section 4). To use this feature,
please set `WithSpikes = FALSE`

as a function parameter in `BASiCS_MCMC`

.

**Important**: BASiCS has been designed in the context of supervised experiments
where the groups of cells (e.g. experimental conditions, cell types) under study
are known a priori (e.g. case-control studies). Therefore, we DO NOT advise the
use of BASiCS in unsupervised settings where the aim is to uncover
sub-populations of cells through clustering.

Parameter estimation is performed using the `BASiCS_MCMC`

function. Downstream
analyses implemented in BASiCS rely on appropriate post-processing of the
output returned by `BASiCS_MCMC`

. Essential parameters for running
`BASiCS_MCMC`

are:

`Data`

: a`SingleCellExperiment`

object created as in Section 3.1`N`

: total number of iterations`Thin`

: length of the thining period (i.e. only every`Thin`

iterations will be stored in the output of the`BASiCS_MCMC`

)`Burn`

: length of burn-in period (i.e. the initial`Burn`

iterations that will be discarded from the output of the`BASiCS_MCMC`

)`Regression`

: if this parameter is set equal to`TRUE`

, the regression BASiCS model will be used (Eling et al. 2017). The latter infers a global regression trend between mean expression and over-dispersion parameters. This trend is subsequently used to derive a*residual over-dispersion*measure that is defined as departures with respect to this trend.**This is now the recommended default setting for BASiCS***.

If the optional parameter `PrintProgress`

is set to `TRUE`

, the R
console will display the progress of the MCMC algorithm.
For other optional parameters refer to `help(BASiCS_MCMC)`

.

Here, we illustrate the usage of `BASiCS_MCMC`

using a built-in
synthetic dataset.

```
Data <- makeExampleBASiCS_Data()
Chain <- BASiCS_MCMC(Data = Data, N = 1000, Thin = 10, Burn = 500,
PrintProgress = FALSE, Regression = TRUE)
```

```
## -----------------------------------------------------
## MCMC sampler has been started: 1000 iterations to go.
## -----------------------------------------------------
## -----------------------------------------------------
## End of Burn-in period.
## -----------------------------------------------------
##
## -----------------------------------------------------
## -----------------------------------------------------
## All 1000 MCMC iterations have been completed.
## -----------------------------------------------------
## -----------------------------------------------------
##
## -----------------------------------------------------
## Please see below a summary of the overall acceptance rates.
## -----------------------------------------------------
##
## Minimum acceptance rate among mu[i]'s: 0.454
## Average acceptance rate among mu[i]'s: 0.66412
## Maximum acceptance rate among mu[i]'s: 0.826
##
##
## Minimum acceptance rate among delta[i]'s: 0.478
## Average acceptance rate among delta[i]'s: 0.55212
## Maximum acceptance rate among delta[i]'s: 0.626
##
##
## Acceptance rate for phi (joint): 0.908
##
##
## Minimum acceptance rate among nu[j]'s: 0.428
## Average acceptance rate among nu[j]'s: 0.5396
## Maximum acceptance rate among nu[j]'s: 0.72
##
##
## Minimum acceptance rate among theta[k]'s: 0.786
## Average acceptance rate among theta[k]'s: 0.786
## Maximum acceptance rate among theta[k]'s: 0.786
##
## -----------------------------------------------------
##
```

As a default, the code above runs the original implementation mode of BASiCS
(spikes without regression; see Section 4).
To use the regression BASiCS model (Eling et al. 2017), please set
`Regression = TRUE`

. To use the no-spikes implementation of BASiCS, please add
`WithSpikes = FALSE`

as an additional parameter.

The `Data`

and `Chain`

(a `BASiCS_Chain`

object) objects created by the code
above can be use for subsequent downstream analyses. See Section
3.2 for highly/lowly variable gene
detection (single group of cells, see also functions `BASiCS_DetectHVG`

and
`BASiCS_DetectLVG`

) and Section 3.3 for
differential expression analyses (two groups of cells, see also function
`BASiCS_TestDE`

).

**Important remarks:**

Please ensure the acceptance rates displayed in the console output of

`BASiCS_MCMC`

are around 0.44. If they are too far from this value, you should increase`N`

and`Burn`

.It is

**essential**to assess the convergence of the MCMC algorithm**before**performing downstream analyses. For guidance regarding this step, refer to the ‘Convergence assessment’ section of this vignette

Typically, setting `N=20000`

, `Thin=20`

and `Burn=10000`

leads to
stable results.

The input dataset for BASiCS must be stored as an `SingleCellExperiment`

object (see *SingleCellExperiment* package).

The `newBASiCS_Data`

function can be used to create the input data object based
on the following information:

`Counts`

: a matrix of raw expression counts with dimensions \(q\) times \(n\). Within this matrix, \(q_0\) rows must correspond to biological genes and \(q-q_0\) rows must correspond to technical spike-in genes. Gene names must be stored as`rownames(Counts)`

.`Tech`

: a logical vector (`TRUE`

/`FALSE`

) with \(q\) elements. If`Tech[i] = FALSE`

the gene`i`

is biological; otherwise the gene is spike-in. This vector must be specified in the same order of genes as in the`Counts`

matrix.`SpikeInfo`

(optional): a`data.frame`

with \(q-q_0\) rows. First column must contain the names associated to the spike-in genes (as in`rownames(Counts)`

). Second column must contain the input number of molecules for the spike-in genes (amount per cell). If a value for this parameter is not provided when calling`newBASiCS_Data`

,`SpikeInfo`

is set as`NULL`

as a default value. In those cases, the`BatchInfo`

argument has to be provided to allow using the no-spikes implementation of BASiCS.`BatchInfo`

(optional): vector of length \(n\) to indicate batch structure (whenever cells have been processed using multiple batches). If a value for this parameter is not provided when calling`newBASiCS_Data`

, BASiCS will assume the data contains a single batch of samples.

For example, the following code generates a synthetic dataset with 50 genes (40 biological and 10 spike-in) and 40 cells.

```
set.seed(1)
Counts <- matrix(rpois(50*40, 2), ncol = 40)
rownames(Counts) <- c(paste0("Gene", 1:40), paste0("Spike", 1:10))
Tech <- c(rep(FALSE,40),rep(TRUE,10))
set.seed(2)
SpikeInput <- rgamma(10,1,1)
SpikeInfo <- data.frame("SpikeID" = paste0("Spike", 1:10),
"SpikeInput" = SpikeInput)
# No batch structure
DataExample <- newBASiCS_Data(Counts, Tech, SpikeInfo)
# With batch structure
DataExample <- newBASiCS_Data(Counts, Tech, SpikeInfo,
BatchInfo = rep(c(1,2), each = 20))
```

Note: BASiCS assumes that a pre-processing quality control step has been applied to remove cells with poor quality data and/or lowly expressed genes that were undetected through sequencing. When analysing multiple groups of cells, the gene filtering step must be jointly applied across all groups to ensure the same genes are retained.

The function `BASiCS_Filter`

can be used to perform this task. For examples,
refer to `help(BASiCS_Filter)`

. Moreover, the *scater* package
provides enhanced functionality for the pre-processing of scRNA-seq datasets.

To convert an existing `SingleCellExperiment`

object (`Data`

) into one that can
be used within BASiCS, meta-information must be stored in the object.

`isSpike(Data, "ERCC") <- Tech`

: the logical vector indicating biological/technical genes (see above) must be stored in the`int_metadata`

slot via the`isSpike`

function.`metadata(Data)`

: the`SpikeInfo`

and`BatchInfo`

objects are stored in the`metadata`

slot of the`SingleCellExperiment`

object:`metadata(Data) <- list(SpikeInput = SpikeInfo[,2], BatchInfo = BatchInfo)`

. Once the additional information is included, the object can be used within BASiCS.

The generation of the input `SingleCellExperiment`

object requires a matrix of
raw counts `Counts`

(columns: cells, rows: genes) after quality control
(e.g. removing low quality cells) and filtering of lowly expressed genes. If
spike-in molecules are contained in `Counts`

, a logical vector `Tech`

is
required to indicate which rows contain technical spike-in molecules and a
`data.frame`

object `SpikeInfo`

containing the names of the spike-in molecules
in the first column and the absolute number of molecules per well in the second
column. More details are provided in section 3.1. If
spike-ins are not available, a vector `BatchInfo`

containing batch information
is required.

We illustrate this analysis using a small extract from the MCMC chain obtained
in (Vallejos, Richardson, and Marioni 2016) when analysing the single cell samples provided in
(Grün, Kester, and Oudenaarden 2014). This is included within `BASiCS`

as the `ChainSC`

dataset.

`data(ChainSC)`

The following code is used to identify **highly variable genes (HVG)** and
**lowly variable genes (LVG)** among these cells. The `VarThreshold`

parameter
sets a lower threshold for the proportion of variability that is assigned to
the biological component (`Sigma`

). In the examples below:

- HVG are defined as those genes for which
**at least**60% of their total variability is attributed to the biological variability component. - LVG are defined as those genes for which
**at most**40% of their total variability is attributed to the biological variability component.

For each gene, these functions return posterior probabilities as a measure of HVG/LVG evidence. A cut-off value for these posterior probabilities is set by controlling the EFDR (as a default option, EFDR is set as 0.10).

```
par(mfrow = c(2,2))
HVG <- BASiCS_DetectHVG(ChainSC, VarThreshold = 0.6, Plot = TRUE)
LVG <- BASiCS_DetectLVG(ChainSC, VarThreshold = 0.2, Plot = TRUE)
```

To access the results of these tests, please use.

`head(HVG$Table)`

```
## GeneIndex GeneName Mu Delta Sigma Prob HVG
## 286 286 Rhox13 9.138169 2.417712 0.7706052 0.9733333 TRUE
## 250 250 Phlda2 9.372955 2.171966 0.7462231 0.9600000 TRUE
## 21 21 Amacr 7.164534 1.598386 0.6652893 0.6933333 TRUE
## 320 320 Smoc1 5.623408 1.805558 0.6700561 0.6800000 FALSE
## 122 122 Engase 7.082273 1.579305 0.6445464 0.6533333 FALSE
## 66 66 Cep170 3.712116 1.580848 0.6261821 0.5600000 FALSE
```

`head(LVG$Table)`

```
## GeneIndex GeneName Mu Delta Sigma Prob LVG
## 20 20 Ahsa1 167.4669 0.06774259 0.09381064 1 TRUE
## 37 37 Atp5g2 344.9518 0.06381609 0.09037440 1 TRUE
## 55 55 Btf3 263.2987 0.04880375 0.06861493 1 TRUE
## 69 69 Chchd2 576.1569 0.03078477 0.04703723 1 TRUE
## 141 141 Fkbp4 334.4190 0.07084192 0.09782010 1 TRUE
## 147 147 Ftl1 2296.2110 0.04504962 0.06656222 1 TRUE
```

**Note**: this decision rule implemented in this function has changed with
respect to the original release of BASiCS (where `EviThreshold`

was defined
such that EFDR = EFNR). However, the new choice is more stable (sometimes, it
was not posible to find a threshold such that EFDR = EFNR).

To illustrate the use of the differential mean expression and differential
over-dispersion tests between two cell populations, we use extracts from the
MCMC chains obtained in (Vallejos, Richardson, and Marioni 2016) when analysing the
(Grün, Kester, and Oudenaarden 2014) dataset (single cells vs pool-and-split samples). These
were obtained by independently running the `BASiCS_MCMC`

function for each
group of cells.

```
data(ChainSC)
data(ChainRNA)
```

```
Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA,
GroupLabel1 = "SC", GroupLabel2 = "PaS",
EpsilonM = log2(1.5), EpsilonD = log2(1.5),
EFDR_M = 0.10, EFDR_D = 0.10,
Offset = TRUE, PlotOffset = FALSE, Plot = TRUE)
```

In `BASiCS_TestDE`

, `EpsilonM`

sets the log2 fold change (log2FC) in expression
(\(\mu\)) and `EpsilonD`

the log2FC in over-dispersion (\(\delta\)). As a default
option: `EpsilonM = EpsilonD = log2(1.5)`

(i.e. the test is set to detect
absolute increases of 50% or above). To adjust for differences in overall mRNA
content, an internal offset correction is performed when `OffSet=TRUE`

.
This is the recommended default setting.

The resulting output list can be displayed using

`head(Test$TableMean)`

```
## GeneName MeanOverall Mean1 Mean2 MeanFC MeanLog2FC ProbDiffMean
## 1 1700094D03Rik 58.539 54.352 62.726 0.854 -0.228 0.000
## 2 1700097N02Rik 37.421 36.919 37.922 0.946 -0.081 0.013
## 3 1810026B05Rik 14.738 10.459 19.018 0.558 -0.842 0.867
## 4 2310008H04Rik 18.844 15.267 22.422 0.672 -0.573 0.493
## 5 2410137M14Rik 19.050 16.674 21.426 0.768 -0.380 0.280
## 6 4930402H24Rik 962.255 967.739 956.771 1.004 0.006 0.000
## ResultDiffMean
## 1 NoDiff
## 2 NoDiff
## 3 PaS+
## 4 NoDiff
## 5 NoDiff
## 6 NoDiff
```

`head(Test$TableDisp)`

```
## GeneName MeanOverall DispOverall Disp1 Disp2 DispFC DispLog2FC
## 1 1700094D03Rik 58.539 0.225 0.321 0.128 2.307 1.206
## 2 1700097N02Rik 37.421 0.338 0.476 0.200 2.342 1.228
## 3 1810026B05Rik 14.738 0.359 0.434 0.284 1.378 0.463
## 4 2310008H04Rik 18.844 0.392 0.441 0.343 1.204 0.268
## 5 2410137M14Rik 19.050 0.396 0.491 0.302 1.663 0.734
## 6 4930402H24Rik 962.255 0.103 0.186 0.020 9.209 3.203
## ProbDiffDisp ResultDiffDisp
## 1 0.880 SC+
## 2 0.813 SC+
## 3 0.520 ExcludedFromTesting
## 4 0.453 NoDiff
## 5 0.640 NoDiff
## 6 1.000 SC+
```

Due to the confounding between mean and over-dispersion that is
typically observed in scRNA-seq datasets, the non-regression BASiCS model
(run using `Regression = FALSE`

as a function parameter in `BASiCS_MCMC`

)
can only be used to assess changes in over-dispersion for those genes in which
the mean expression does not change between the groups. In this case, we
recommend users to use `EpsilonM = 0`

as a conservative option to avoid
changes in over-dispersion to be confounded by mean expression (the genes for
which mean expression changes are marked as `ExcludedFromTesting`

in the
`Test$TableDisp$ResultDiffDisp`

slot).

```
Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA,
GroupLabel1 = "SC", GroupLabel2 = "PaS",
EpsilonM = 0, EpsilonD = log2(1.5),
EFDR_M = 0.10, EFDR_D = 0.10,
Offset = TRUE, PlotOffset = FALSE, Plot = FALSE)
```

**Note**: If the regression BASiCS model has been
used (`Regression = TRUE`

as a function parameter in `BASiCS_MCMC`

),
`BASiCS_TestDE`

will also report changes in residual over-dispersion
(not confounded by mean expression) between the groups (see Section
4 in this vignette).

Beyond its original implementation, BASiCS has been extended to adress experimental designs in which spike-in molecules are not available as well as to address the confounding that is typically observed between mean and over-dispersion for scRNA-seq datasets (Eling et al. 2017). Alternative implementation modes are summarised below:

As a default, the `BASiCS_MCMC`

function uses `WithSpikes = TRUE`

.

`WithSpikes = FALSE`

When technical spike-in genes are not available, BASiCS uses a horizontal
integration strategy which borrows information across multiple technical
replicates (Eling et al. 2017). Therefore, `BASiCS_MCMC`

will fail to run if a
single batch of samples is provided. **Note:** batch information must be
provided via the `BatchInfo`

argument when using the `newBASiCS_Data`

function.
Additionally, when creating the `BASiCS_Data`

object, the `SpikeInfo`

slot must
be set as `SpikeInfo = NULL`

.

```
DataNoSpikes <- newBASiCS_Data(Counts, Tech, SpikeInfo = NULL,
BatchInfo = rep(c(1,2), each = 20))
ChainNoSpikes <- BASiCS_MCMC(Data = DataNoSpikes, N = 1000,
Thin = 10, Burn = 500,
WithSpikes = FALSE, Regression = TRUE,
PrintProgress = FALSE)
```

```
## -----------------------------------------------------
## MCMC sampler has been started: 1000 iterations to go.
## -----------------------------------------------------
## -----------------------------------------------------
## End of Burn-in period.
## -----------------------------------------------------
##
## -----------------------------------------------------
## -----------------------------------------------------
## All 1000 MCMC iterations have been completed.
## -----------------------------------------------------
## -----------------------------------------------------
##
## -----------------------------------------------------
## Please see below a summary of the overall acceptance rates.
## -----------------------------------------------------
##
## Minimum acceptance rate among mu[i]'s: 0.332
## Average acceptance rate among mu[i]'s: 0.47664
## Maximum acceptance rate among mu[i]'s: 0.67
##
##
## Minimum acceptance rate among delta[i]'s: 0.706
## Average acceptance rate among delta[i]'s: 0.7634
## Maximum acceptance rate among delta[i]'s: 0.822
##
##
## Minimum acceptance rate among nu[jk]'s: 0.912
## Average acceptance rate among nu[jk]'s: 0.9589
## Maximum acceptance rate among nu[jk]'s: 0.984
##
##
## Minimum acceptance rate among theta[k]'s: 0.79
## Average acceptance rate among theta[k]'s: 0.793
## Maximum acceptance rate among theta[k]'s: 0.796
##
##
## -----------------------------------------------------
##
```

`Regression = TRUE`

The BASiCS model uses a joint informative prior formulation to account for the
relationship between mean and over-dispersion gene-specific parameters. The
latter is used to infer a global regression trend between these parameters and,
subsequently, to derive a *residual over-dispersion* measure that is
defined as departures with respect to this trend.

```
DataRegression <- newBASiCS_Data(Counts, Tech, SpikeInfo,
BatchInfo = rep(c(1,2), each = 20))
ChainRegression <- BASiCS_MCMC(Data = DataRegression, N = 1000,
Thin = 10, Burn = 500,
Regression = TRUE,
PrintProgress = FALSE)
```

```
## -----------------------------------------------------
## MCMC sampler has been started: 1000 iterations to go.
## -----------------------------------------------------
## -----------------------------------------------------
## End of Burn-in period.
## -----------------------------------------------------
##
## -----------------------------------------------------
## -----------------------------------------------------
## All 1000 MCMC iterations have been completed.
## -----------------------------------------------------
## -----------------------------------------------------
##
## -----------------------------------------------------
## Please see below a summary of the overall acceptance rates.
## -----------------------------------------------------
##
## Minimum acceptance rate among mu[i]'s: 0.406
## Average acceptance rate among mu[i]'s: 0.498
## Maximum acceptance rate among mu[i]'s: 0.602
##
##
## Minimum acceptance rate among delta[i]'s: 0.74
## Average acceptance rate among delta[i]'s: 0.79905
## Maximum acceptance rate among delta[i]'s: 0.838
##
##
## Acceptance rate for phi (joint): 0.84
##
##
## Minimum acceptance rate among nu[j]'s: 0.776
## Average acceptance rate among nu[j]'s: 0.84355
## Maximum acceptance rate among nu[j]'s: 0.958
##
##
## Minimum acceptance rate among theta[k]'s: 0.734
## Average acceptance rate among theta[k]'s: 0.735
## Maximum acceptance rate among theta[k]'s: 0.736
##
## -----------------------------------------------------
##
```

This implementation provides additional functionality when performing
downstream analyses. These are illustrated below using a small extract from
the MCMC chain obtained when analysing the dataset provided in
(Grün, Kester, and Oudenaarden 2014) (single cell versus pool-and-split samples). These are
included within `BASiCS`

as the `ChainSCReg`

and `ChainRNAReg`

datasets.

To visualize the fit between over-dispersion \(\delta_i\) and mean expression $ _i$ the following function can be used.

```
data("ChainRNAReg")
BASiCS_showFit(ChainRNAReg)
```

The `BASiCS_TestDE`

test function will automatically perform differential
variability testing based on the residual over-dispersion parameters
\(\epsilon_i\) when its output includes two `Chain`

objects that were generated
by the regression BASiCS model.

```
data("ChainSCReg")
Test <- BASiCS_TestDE(Chain1 = ChainSCReg, Chain2 = ChainRNAReg,
GroupLabel1 = "SC", GroupLabel2 = "PaS",
EpsilonM = log2(1.5), EpsilonD = log2(1.5),
EpsilonR = log2(1.5)/log2(exp(1)),
EFDR_M = 0.10, EFDR_D = 0.10,
Offset = TRUE, PlotOffset = FALSE, Plot = FALSE)
```

This test function outputs an extra slot containing the results of the
differential testing residual over-dispersion test. Only genes that are
expressed in at least 2 cells (in both groups) are included in the test.
Genes that do not satisfy this condition are marked as `ExcludedFromRegression`

in the `Test$TableResDisp$ResultDiffResDisp`

slot. By performing the regression,
all genes can be tested for changes in expression variability independent of
changes in mean expression.

`head(Test$TableResDisp, n = 2)`

```
## GeneName MeanOverall ResDispOverall ResDisp1 ResDisp2 ResDispDistance
## 1 1700094D03Rik 57.353 0.100 0.420 -0.221 0.632
## 2 1700097N02Rik 36.971 0.152 0.438 -0.133 0.432
## ProbDiffResDisp ResultDiffResDisp
## 1 0.707 NoDiff
## 2 0.613 NoDiff
```

**Note:** Additional parameters for this sampler include: `k`

number of
regression components (`k`

-2 radial basis functions, one intercept and one
linear component), `Var`

the scale parameter influencing the width of the basis
functions and `eta`

the degrees of freedom. For additional details about these
choices, please refer to Eling et al. (2017).

To externally store the output of `BASiCS_MCMC`

(recommended), additional
parameters `StoreChains`

, `StoreDir`

and `RunName`

are required. For example:

```
Data <- makeExampleBASiCS_Data()
Chain <- BASiCS_MCMC(Data, N = 1000, Thin = 10, Burn = 500, Regression = TRUE,
PrintProgress = FALSE, StoreChains = TRUE,
StoreDir = tempdir(), RunName = "Example")
```

```
## -----------------------------------------------------
## MCMC sampler has been started: 1000 iterations to go.
## -----------------------------------------------------
## -----------------------------------------------------
## End of Burn-in period.
## -----------------------------------------------------
##
## -----------------------------------------------------
## -----------------------------------------------------
## All 1000 MCMC iterations have been completed.
## -----------------------------------------------------
## -----------------------------------------------------
##
## -----------------------------------------------------
## Please see below a summary of the overall acceptance rates.
## -----------------------------------------------------
##
## Minimum acceptance rate among mu[i]'s: 0.454
## Average acceptance rate among mu[i]'s: 0.66412
## Maximum acceptance rate among mu[i]'s: 0.826
##
##
## Minimum acceptance rate among delta[i]'s: 0.478
## Average acceptance rate among delta[i]'s: 0.55212
## Maximum acceptance rate among delta[i]'s: 0.626
##
##
## Acceptance rate for phi (joint): 0.908
##
##
## Minimum acceptance rate among nu[j]'s: 0.428
## Average acceptance rate among nu[j]'s: 0.5396
## Maximum acceptance rate among nu[j]'s: 0.72
##
##
## Minimum acceptance rate among theta[k]'s: 0.786
## Average acceptance rate among theta[k]'s: 0.786
## Maximum acceptance rate among theta[k]'s: 0.786
##
## -----------------------------------------------------
##
```

In this example, the output of `BASiCS_MCMC`

will be stored as a `BASiCS_Chain`

object in the file “chain_Example.Rds”, within the `tempdir()`

directory.

To load pre-computed MCMC chains,

`Chain <- BASiCS_LoadChain("Example", StoreDir = tempdir()) `

To assess convergence of the chain, the convergence diagnostics provided by the
package `coda`

can be used. Additionally, the chains can be visually inspected.
For example:

`plot(Chain, Param = "mu", Gene = 1, log = "y")`