library(scDesign3)
library(SingleCellExperiment)
library(ggplot2)
theme_set(theme_bw())

0.1 Introduction

scDesign3 is a unified probabilistic framework that generates realistic in silico high-dimensional single-cell omics data of various cell states, including discrete cell types, continuous trajectories, and spatial locations by learning from real datasets. Since the functions of scDesign3 is very comprehensive, here we only introduce how scDesign3 simulates an scRNA-seq dataset with one continuous developmental trajectory. For more information, please check the Articles on our website: (https://songdongyuan1994.github.io/scDesign3/docs/index.html).

0.2 Read in the reference data

The raw data is from the scvelo, which describes pancreatic endocrinogenesis. We pre-select the top 1000 highly variable genes and filter out some cell types to ensure a single trajectory.

example_sce <- readRDS((url("https://figshare.com/ndownloader/files/40581992")))
print(example_sce)
#> class: SingleCellExperiment 
#> dim: 1000 2087 
#> metadata(5): clusters_coarse_colors clusters_colors day_colors
#>   neighbors pca
#> assays(6): X spliced ... cpm logcounts
#> rownames(1000): Pyy Iapp ... Eya2 Kif21a
#> rowData names(1): highly_variable_genes
#> colnames(2087): AAACCTGAGAGGGATA AAACCTGGTAAGTGGC ... TTTGTCAAGTGACATA
#>   TTTGTCAAGTGTGGCA
#> colData names(7): clusters_coarse clusters ... sizeFactor pseudotime
#> reducedDimNames(4): X_pca X_umap PCA UMAP
#> mainExpName: NULL
#> altExpNames(0):

To save computational time, we only use the top 100 genes.

example_sce <- example_sce[1:100, ]

0.3 Simulation

The function scdesign3() takes in a SinglecellExperiment object with the cell covariates (such as cell types, pesudotime, or spatial coordinates) stored in the colData of the SinglecellExperiment object.

set.seed(123)
example_simu <- scdesign3(
    sce = example_sce,
    assay_use = "counts",
    celltype = "cell_type",
    pseudotime = "pseudotime",
    spatial = NULL,
    other_covariates = NULL,
    mu_formula = "s(pseudotime, k = 10, bs = 'cr')",
    sigma_formula = "1", # If you want your dispersion also varies along pseudotime, use "s(pseudotime, k = 5, bs = 'cr')"
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    DT = TRUE,
    pseudo_obs = FALSE,
    return_model = FALSE,
    nonzerovar = FALSE
  )

The output of scdesign3() is a list which includes:

  • new_count: This is the synthetic count matrix generated by scdesign3().
  • new_covariate:
    • If the parameter ncell is set to a number that is different from the number of cells in the input data, this will be a matrix that has the new cell covariates that are used for generating new data.
    • If the parameter ncell is the default value, this will be NULL.
  • model_aic: This is a vector include the genes’ marginal models’ AIC, fitted copula’s AIC, and total AIC, which is the sum of the previous two AIC.
  • model_bic: This is a vector include the genes’ marginal models’ BIC, fitted copula’s BIC, and total BIC, which is the sum of the previous two BIC.
  • marginal_list:
    • If the parameter return_model is set to TRUE, this will be a list which contains the fitted gam or gamlss model for all genes in the input data.
    • If the parameter return_model is set to the default value FALSE, this will be NULL.
  • corr_list:
    • If the parameter return_model is set to TRUE, this will be a list which contains the either a correlation matrix (when copula = "gaussian") or the fitted vine copula (when copula = "vine) for each user specified correlation groups (based on the parameter corr_by).
    • If the parameter return_model is set to the default value FALSE, this will be NULL.

In this example, since we did not change the parameter ncell, the synthetic count matrix will have the same dimension as the input count matrix.

dim(example_simu$new_count)
#> [1]  100 2087

Then, we can create the SinglecellExperiment object using the synthetic count matrix and store the logcounts to the input and synthetic SinglecellExperiment objects.

logcounts(example_sce) <- log1p(counts(example_sce))
simu_sce <- SingleCellExperiment(list(counts = example_simu$new_count), colData = example_simu$new_covariate)
logcounts(simu_sce) <- log1p(counts(simu_sce))

0.4 Visualization

set.seed(123)
compare_figure <- plot_reduceddim(ref_sce = example_sce, 
                                  sce_list = list(simu_sce), 
                                  name_vec = c("Reference", "scDesign3"),
                                  assay_use = "logcounts", 
                                  if_plot = TRUE, 
                                  color_by = "pseudotime", 
                                  n_pc = 20)
plot(compare_figure$p_umap)

0.5 Session information

sessionInfo()
#> R version 4.3.2 Patched (2023-11-13 r85521)
#> 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] ggplot2_3.4.4               SingleCellExperiment_1.24.0
#>  [3] SummarizedExperiment_1.32.0 Biobase_2.62.0             
#>  [5] GenomicRanges_1.54.1        GenomeInfoDb_1.38.5        
#>  [7] IRanges_2.36.0              S4Vectors_0.40.2           
#>  [9] BiocGenerics_0.48.1         MatrixGenerics_1.14.0      
#> [11] matrixStats_1.2.0           scDesign3_1.0.1            
#> [13] BiocStyle_2.30.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.4            xfun_0.41               bslib_0.6.1            
#>  [4] lattice_0.22-5          gamlss_5.4-20           vctrs_0.6.5            
#>  [7] tools_4.3.2             bitops_1.0-7            generics_0.1.3         
#> [10] parallel_4.3.2          tibble_3.2.1            fansi_1.0.6            
#> [13] highr_0.10              pkgconfig_2.0.3         Matrix_1.6-5           
#> [16] lifecycle_1.0.4         GenomeInfoDbData_1.2.11 farver_2.1.1           
#> [19] compiler_4.3.2          munsell_0.5.0           htmltools_0.5.7        
#> [22] sass_0.4.8              RCurl_1.98-1.14         yaml_2.3.8             
#> [25] pillar_1.9.0            crayon_1.5.2            jquerylib_0.1.4        
#> [28] MASS_7.3-60.0.1         openssl_2.1.1           cachem_1.0.8           
#> [31] DelayedArray_0.28.0     viridis_0.6.4           abind_1.4-5            
#> [34] mclust_6.0.1            RSpectra_0.16-1         nlme_3.1-164           
#> [37] tidyselect_1.2.0        digest_0.6.34           mvtnorm_1.2-4          
#> [40] dplyr_1.1.4             bookdown_0.37           labeling_0.4.3         
#> [43] splines_4.3.2           gamlss.dist_6.1-1       fastmap_1.1.1          
#> [46] grid_4.3.2              gamlss.data_6.0-2       colorspace_2.1-0       
#> [49] cli_3.6.2               SparseArray_1.2.3       magrittr_2.0.3         
#> [52] S4Arrays_1.2.0          survival_3.5-7          utf8_1.2.4             
#> [55] withr_2.5.2             scales_1.3.0            rmarkdown_2.25         
#> [58] XVector_0.42.0          umap_0.2.10.0           gridExtra_2.3          
#> [61] reticulate_1.34.0       png_0.1-8               askpass_1.2.0          
#> [64] evaluate_0.23           knitr_1.45              viridisLite_0.4.2      
#> [67] irlba_2.3.5.1           mgcv_1.9-1              rlang_1.1.3            
#> [70] Rcpp_1.0.12             glue_1.7.0              BiocManager_1.30.22    
#> [73] jsonlite_1.8.8          R6_2.5.1                zlibbioc_1.48.0