An introduction to BASiCStan

A Stan implementation of BASiCS

Outline

This package implements alternative inference methods for BASiCS. The original package uses adaptive Metropolis within Gibbs sampling, while BASiCStan uses Stan to enable the use of maximum a posteriori estimation, variational inference, and Hamiltonian Monte Carlo. These each have advantages for different use cases.

Use

To use BASiCStan, we can first use BASICS_MockSCE() to generate an toy example dataset. We will also set a seed for reproducibility.

suppressPackageStartupMessages(library("BASiCStan"))
set.seed(42)
sce <- BASiCS_MockSCE()

The interface for running MCMC using BASiCS and using alternative inference methods using Stan is very similar. It is worth noting that the joint prior between mean and over-dispersion parameters, corresponding to Regression = TRUE, is the only supported mode in BASiCStan().

amwg_fit <- BASiCS_MCMC(
    sce,
    N = 200,
    Thin = 10,
    Burn = 100,
    Regression = TRUE
)
#> altExp 'spike-ins' is assumed to contain spike-in genes.
#> see help(altExp) for details.
#> Running with spikes BASiCS sampler (regression case) ...
#> -------------------------------------------------------------
#> MCMC running time 
#> -------------------------------------------------------------
#> user: 0.343
#> system: 0
#> elapsed: 0.343
#> -------------------------------------------------------------
#> Output 
#> -------------------------------------------------------------
#> -------------------------------------------------------------
#> BASiCS version 2.14.0 : 
#> vertical integration (spikes case) 
#> -------------------------------------------------------------
stan_fit <- BASiCStan(sce, Method = "sampling", iter = 10)
#> Warning: There were 5 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems

The output of BASiCStan() is an object of class BASiCS_Chain, similar to the output of BASiCS_MCMC(). Therefore, you could use these as you would an object generated using Metropolis within Gibbs sampling. For example, we can plot the relationship between mean and over-dispersion estimated using the joint regression prior:

BASiCS_ShowFit(amwg_fit)

BASiCS_ShowFit(stan_fit)

Stan MCMC diagnostics

Using Stan has advantages outside of the variety of inference engines available. By returning a Stan object that we can later convert to a BASiCS_Chain object, we can leverage an even broader set of functionality. For example, Stan has the ability to easily generate MCMC diagnostics using simple functions. For example, summary() outputs a number of useful per-chain and across-chain diagnostics:

stan_obj <- BASiCStan(sce, ReturnBASiCS = FALSE, Method = "sampling", iter = 10)
#> 
#> SAMPLING FOR MODEL 'basics_regression' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.003532 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 35.32 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: No variance estimation is
#> Chain 1:          performed for num_warmup < 20
#> Chain 1: 
#> Chain 1: Iteration: 1 / 10 [ 10%]  (Warmup)
#> Chain 1: Iteration: 2 / 10 [ 20%]  (Warmup)
#> Chain 1: Iteration: 3 / 10 [ 30%]  (Warmup)
#> Chain 1: Iteration: 4 / 10 [ 40%]  (Warmup)
#> Chain 1: Iteration: 5 / 10 [ 50%]  (Warmup)
#> Chain 1: Iteration: 6 / 10 [ 60%]  (Sampling)
#> Chain 1: Iteration: 7 / 10 [ 70%]  (Sampling)
#> Chain 1: Iteration: 8 / 10 [ 80%]  (Sampling)
#> Chain 1: Iteration: 9 / 10 [ 90%]  (Sampling)
#> Chain 1: Iteration: 10 / 10 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.054 seconds (Warm-up)
#> Chain 1:                0.034 seconds (Sampling)
#> Chain 1:                0.088 seconds (Total)
#> Chain 1:
#> Warning: There were 5 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
head(summary(stan_obj)$summary)
#>               mean se_mean sd     2.5%      25%      50%      75%    97.5%
#> log_mu[1] 3.098420     NaN  0 3.098420 3.098420 3.098420 3.098420 3.098420
#> log_mu[2] 1.058488     NaN  0 1.058488 1.058488 1.058488 1.058488 1.058488
#> log_mu[3] 2.107672     NaN  0 2.107672 2.107672 2.107672 2.107672 2.107672
#> log_mu[4] 2.245089     NaN  0 2.245089 2.245089 2.245089 2.245089 2.245089
#> log_mu[5] 2.002693     NaN  0 2.002693 2.002693 2.002693 2.002693 2.002693
#> log_mu[6] 1.457520     NaN  0 1.457520 1.457520 1.457520 1.457520 1.457520
#>           n_eff Rhat
#> log_mu[1]   NaN  NaN
#> log_mu[2]   NaN  NaN
#> log_mu[3]   NaN  NaN
#> log_mu[4]   NaN  NaN
#> log_mu[5]   NaN  NaN
#> log_mu[6]   NaN  NaN

We can then convert this object to a BASiCS_Chain and carry on a workflow as before:

Stan2BASiCS(stan_obj)
#> An object of class BASiCS_Chain
#>  5 MCMC samples.
#>  Dataset contains 100 biological genes and 100 cells (2 batches). 
#>  Object stored using BASiCS version:  2.14.0 
#>  Parameters:  mu delta s nu theta epsilon phi beta

Session info

sessionInfo()
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] BASiCStan_1.4.0             rstan_2.32.3               
#>  [3] StanHeaders_2.26.28         BASiCS_2.14.0              
#>  [5] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
#>  [7] Biobase_2.62.0              GenomicRanges_1.54.0       
#>  [9] GenomeInfoDb_1.38.0         IRanges_2.36.0             
#> [11] S4Vectors_0.40.0            BiocGenerics_0.48.0        
#> [13] MatrixGenerics_1.14.0       matrixStats_1.0.0          
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.0-7              gridExtra_2.3            
#>   [3] inline_0.3.19             rlang_1.1.1              
#>   [5] magrittr_2.0.3            compiler_4.3.1           
#>   [7] loo_2.6.0                 DelayedMatrixStats_1.24.0
#>   [9] callr_3.7.3               vctrs_0.6.4              
#>  [11] pkgconfig_2.0.3           crayon_1.5.2             
#>  [13] fastmap_1.1.1             backports_1.4.1          
#>  [15] XVector_0.42.0            ellipsis_0.3.2           
#>  [17] labeling_0.4.3            scuttle_1.12.0           
#>  [19] utf8_1.2.4                promises_1.2.1           
#>  [21] rmarkdown_2.25            ps_1.7.5                 
#>  [23] xfun_0.40                 bluster_1.12.0           
#>  [25] zlibbioc_1.48.0           cachem_1.0.8             
#>  [27] beachmat_2.18.0           jsonlite_1.8.7           
#>  [29] later_1.3.1               DelayedArray_0.28.0      
#>  [31] BiocParallel_1.36.0       prettyunits_1.2.0        
#>  [33] irlba_2.3.5.1             parallel_4.3.1           
#>  [35] cluster_2.1.4             R6_2.5.1                 
#>  [37] bslib_0.5.1               limma_3.58.0             
#>  [39] jquerylib_0.1.4           Rcpp_1.0.11              
#>  [41] assertthat_0.2.1          knitr_1.44               
#>  [43] splines_4.3.1             httpuv_1.6.12            
#>  [45] Matrix_1.6-1.1            igraph_1.5.1             
#>  [47] tidyselect_1.2.0          abind_1.4-5              
#>  [49] yaml_2.3.7                viridis_0.6.4            
#>  [51] codetools_0.2-19          miniUI_0.1.1.1           
#>  [53] curl_5.1.0                processx_3.8.2           
#>  [55] pkgbuild_1.4.2            lattice_0.22-5           
#>  [57] tibble_3.2.1              withr_2.5.1              
#>  [59] shiny_1.7.5.1             posterior_1.4.1          
#>  [61] coda_0.19-4               evaluate_0.22            
#>  [63] RcppParallel_5.1.7        pillar_1.9.0             
#>  [65] tensorA_0.36.2            checkmate_2.2.0          
#>  [67] distributional_0.3.2      generics_0.1.3           
#>  [69] RCurl_1.98-1.12           ggplot2_3.4.4            
#>  [71] rstantools_2.3.1.1        sparseMatrixStats_1.14.0 
#>  [73] munsell_0.5.0             scales_1.2.1             
#>  [75] xtable_1.8-4              glue_1.6.2               
#>  [77] metapod_1.10.0            tools_4.3.1              
#>  [79] hexbin_1.28.3             BiocNeighbors_1.20.0     
#>  [81] ScaledMatrix_1.10.0       locfit_1.5-9.8           
#>  [83] scran_1.30.0              cowplot_1.1.1            
#>  [85] grid_4.3.1                QuickJSR_1.0.7           
#>  [87] edgeR_4.0.0               colorspace_2.1-0         
#>  [89] GenomeInfoDbData_1.2.11   BiocSingular_1.18.0      
#>  [91] cli_3.6.1                 rsvd_1.0.5               
#>  [93] fansi_1.0.5               viridisLite_0.4.2        
#>  [95] S4Arrays_1.2.0            dplyr_1.1.3              
#>  [97] V8_4.4.0                  glmGamPoi_1.14.0         
#>  [99] gtable_0.3.4              sass_0.4.7               
#> [101] digest_0.6.33             SparseArray_1.2.0        
#> [103] dqrng_0.3.1               farver_2.1.1             
#> [105] htmltools_0.5.6.1         lifecycle_1.0.3          
#> [107] statmod_1.5.0             mime_0.12                
#> [109] ggExtra_0.10.1            MASS_7.3-60