BASiCS (version 0.3.1)

BASiCS_MCMC: 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, ...)

Arguments

Data

A BASiCS_Data-class object.

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.

...

Optional parameters.

PriorParam

List of 7 elements, containing the hyper-parameter values required for the adopted prior (see Vallejos et al, 2015). All elements must be positive real numbers.

s2.mu

Scale hyper-parameter for the log-Normal(0,s2.mu) prior that is shared by all gene-specific expression rate parameters \(\mu[i]\). Default: s2.mu = 1.

a.delta

Shape hyper-parameter for the Gamma(a.delta,b.delta) prior that is shared by all gene-specific biological cell-to-cell heterogeneity hyper-parameters \(\delta[i]\). Default: a.delta = 1.

b.delta

Rate hyper-parameter for the Gamma(a.delta,b.delta) prior that is shared by all gene-specific biological cell-to-cell heterogeneity hyper-parameters \(\delta[i]\). Default: b.delta = 1.

p.phi

Dirichlet hyper-parameter for the joint of all (scaled by n) cell-specific mRNA content normalising constants \(\phi[j] / n\). Default: p.phi = rep(1, n).

a.s

Shape hyper-parameter for the Gamma(a.s,b.s) prior that is shared by all cell-specific capture efficiency normalising constants \(s[j]\). Default: a.s = 1.

b.s

Rate hyper-parameter for the Gamma(a.s,b.s) prior that is shared by all cell-specific capture efficiency normalising constants \(s[j]\). Default: b.s = 1.

a.theta

Shape hyper-parameter for the Gamma(a.theta,b.theta) prior for technical noise hyper-parameter \(\theta\). Default: a.theta = 1.

b.theta

Rate hyper-parameter for the Gamma(a.theta,b.theta) prior for technical noise hyper-parameter \(\theta\). Default: b.theta = 1.

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 = T, the slots of the generated BASiCS_Chain object are stored in separate .txt files. Each row of the output file containing an interation (RunName argument used to index file names). Default: StoreChains = F.

StoreAdapt

If StoreAdapt = T, trajectory of adaptive proposal variances (in log-scale) for each parameter are stored in separate .txt files. Each row of the output file containing an interation (RunName argument used to index file names). Default: StoreAdapt = F.

StoreDir

Directory where output files are stored. Only required if StoreChains = TRUE and/or StoreAdapt = TRUE). Default: StoreDir = getwd().

RunName

String used to index `.txt` files storing chains and/or adaptive proposal variances.

PrintProgress

If PrintProgress = FALSE, console-based progress report is suppressed.

ls.phi0

Starting value for the adaptive concentration parameter of the Metropolis proposals for phi.

Value

An object of class BASiCS_Chain-class.

References

Vallejos, Marioni and Richardson (2015). Bayesian Analysis of Single-Cell Sequencing data.

Examples

# NOT RUN {
# Built-in simulated dataset
Data = makeExampleBASiCS_Data()

# Only a short run of the MCMC algorithm for illustration purposes
# Longer runs migth be required to reach convergence
MCMC_Output <- BASiCS_MCMC(Data, N = 10000, Thin = 10, Burn = 5000, PrintProgress = FALSE)
head(displayChainBASiCS(MCMC_Output, Param = "mu"))
head(displayChainBASiCS(MCMC_Output, Param = "delta"))
head(displayChainBASiCS(MCMC_Output, Param = "phi"))
head(displayChainBASiCS(MCMC_Output, Param = "s"))
head(displayChainBASiCS(MCMC_Output, Param = "nu"))
head(displayChainBASiCS(MCMC_Output, Param = "theta"))

# Traceplots
plot(MCMC_Output, Param = "mu", Gene = 1)
plot(MCMC_Output, Param = "delta", Gene = 1)
plot(MCMC_Output, Param = "phi", Cell = 1)
plot(MCMC_Output, Param = "s", Cell = 1)
plot(MCMC_Output, Param = "nu", Cell = 1)
plot(MCMC_Output, Param = "theta", Batch = 1)

# Calculating posterior medians and 95% HPD intervals
MCMC_Summary <- Summary(MCMC_Output)
head(displaySummaryBASiCS(MCMC_Summary, Param = "mu"))
head(displaySummaryBASiCS(MCMC_Summary, Param = "delta"))
head(displaySummaryBASiCS(MCMC_Summary, Param = "phi"))
head(displaySummaryBASiCS(MCMC_Summary, Param = "s"))
head(displaySummaryBASiCS(MCMC_Summary, Param = "nu"))
head(displaySummaryBASiCS(MCMC_Summary, Param = "theta"))

# Graphical display of posterior medians and 95% HPD intervals
plot(MCMC_Summary, Param = "mu", main = "All genes")
plot(MCMC_Summary, Param = "mu", Genes = 1:10, main = "First 10 genes")
plot(MCMC_Summary, Param = "delta", main = "All genes")
plot(MCMC_Summary, Param = "delta", Genes = c(2,5,10,50,100), main = "5 customized genes")
plot(MCMC_Summary, Param = "phi", main = "All cells")
plot(MCMC_Summary, Param = "phi", Cells = 1:5, main = "First 5 cells")
plot(MCMC_Summary, Param = "s", main = "All cells")
plot(MCMC_Summary, Param = "s", Cells = 1:5, main = "First 5 cells")
plot(MCMC_Summary, Param = "nu", main = "All cells")
plot(MCMC_Summary, Param = "nu", Cells = 1:5, main = "First 5 cells")
plot(MCMC_Summary, Param = "theta")

# To constrast posterior medians of cell-specific parameters
plot(MCMC_Summary, Param = "phi", Param2 = "s")
plot(MCMC_Summary, Param = "phi", Param2 = "nu")
plot(MCMC_Summary, Param = "s", Param2 = "nu")

# To constrast posterior medians of gene-specific parameters
plot(MCMC_Summary, Param = "mu", Param2 = "delta", log = "x")

# Highly and lowly variable genes detection
#DetectHVG <- BASiCS_DetectHVG(Data, MCMC_Output, VarThreshold = 0.70, Plot = TRUE)
#DetectLVG <- BASiCS_DetectLVG(Data, MCMC_Output, VarThreshold = 0.40, Plot = TRUE)

#plot(MCMC_Summary, Param = "mu", Param2 = "delta", log = "x", col = 8)
#points(DetectHVG$Table[DetectHVG$Table$HVG==1,2], DetectHVG$Table[DetectHVG$Table$HVG==1,3],
#       pch = 16, col = "red", cex = 1)
#points(DetectLVG$Table[DetectLVG$Table$LVG==1,2], DetectLVG$Table[DetectLVG$Table$LVG==1,3],
#       pch = 16, col = "blue", cex = 1)

# If variance thresholds are not fixed
#BASiCS_VarThresholdSearchHVG(Data, MCMC_Output, VarThresholdsGrid = seq(0.70,0.75,by=0.01))
#BASiCS_VarThresholdSearchLVG(Data, MCMC_Output, VarThresholdsGrid = seq(0.40,0.45,by=0.01))
# }