scDesign3

https://songdongyuan1994.github.io/scDesign3/docs/
Admin
scDesign3 is an all-in-one single-cell and spatial omics data simulation tool that generates...
Tokens:63,765
Snippets:389
Trust Score:9.7
Update:1 month ago
Show doc for...
Tokens:
Raw
### simu_new example workflow

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/simu_new.html

Complete example of using simu_new with example_sce data. It constructs data, fits marginal and copula models, extracts parameters, and simulates new counts. The output includes progress messages from fit_copula and simu_new.

```R
data(example_sce)
  my_data <- construct_data(
  sce = example_sce,
  assay_use = "counts",
  celltype = "cell_type",
  pseudotime = "pseudotime",
  spatial = NULL,
  other_covariates = NULL,
  corr_by = "1"
  )
  my_marginal <- fit_marginal(
  data = my_data,
  mu_formula = "s(pseudotime, bs = 'cr', k = 10)",
  sigma_formula = "1",
  family_use = "nb",
  n_cores = 1,
  usebam = FALSE
  )
  my_copula <- fit_copula(
  sce = example_sce,
  assay_use = "counts",
  marginal_list = my_marginal,
  family_use = c(rep("nb", 5), rep("zip", 5)),
  copula = "vine",
  n_cores = 1,
  input_data = my_data$dat
  )
#> Convert Residuals to Uniform
#> Converting End
#> Copula group 1 starts
#> Vine Copula Estimation Starts
#> Time difference of 0.1432807 secs
#> Vine Copula Estimation Ends
  my_para <- extract_para(
    sce = example_sce,
    marginal_list = my_marginal,
    n_cores = 1,
    family_use = c(rep("nb", 5), rep("zip", 5)),
    new_covariate = my_data$new_covariate,
    data = my_data$dat
  )
  my_newcount <- simu_new(
  sce = example_sce,
  mean_mat = my_para$mean_mat,
  sigma_mat = my_para$sigma_mat,
  zero_mat = my_para$zero_mat,
  quantile_mat = NULL,
  copula_list = my_copula$copula_list,
  n_cores = 1,
  family_use = c(rep("nb", 5), rep("zip", 5)),
  input_data = my_data$dat,
  new_covariate = my_data$new_covariate,
  important_feature = my_copula$important_feature,
  filtered_gene = my_data$filtered_gene
  )
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
```

--------------------------------

### Example for ba() printing "No example"

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/ba.html

Shows the only example provided for the ba() function, which simply prints the string "No example". The output is displayed as a comment with the result.

```r
print("No example")
#> [1] "No example"
```

--------------------------------

### Construct input data with example_sce

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/construct_data.html

Example usage of construct_data with example_sce, specifying assay_use, celltype, pseudotime, and corr_by. The example does not set ncell, so newCovariate will be NULL.

```r
data(example_sce)
my_data <- construct_data(
  sce = example_sce,
  assay_use = "counts",
  celltype = "cell_type",
  pseudotime = "pseudotime",
  spatial = NULL,
  other_covariates = NULL,
  corr_by = "1"
)
```

--------------------------------

### Run scDesign3 simulation with scdesign3()

Source: https://songdongyuan1994.github.io/scDesign3/docs/index.html

Quick start example of running the scDesign3 simulator. The scdesign3() function takes a SingleCellExperiment object with cell covariates stored in colData and returns a simulated dataset. This example uses a negative binomial family, pseudotime as a smooth term, and a Gaussian copula.

```r
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 = "s(pseudotime, k = 5, bs = 'cr')",
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    fastmvn = FALSE,
    DT = TRUE,
    pseudo_obs = FALSE,
    family_set = c("gauss", "indep"),
    important_feature = "all",
    nonnegative = TRUE,
    return_model = FALSE,
    nonzerovar = FALSE,
    parallelization = "mcmapply",
    BPPARAM = NULL,
    trace = FALSE
  )
```

--------------------------------

### Install scDesign3 from GitHub with devtools

Source: https://songdongyuan1994.github.io/scDesign3/docs/index.html

Installs the development version of scDesign3 from GitHub using devtools. Requires the devtools package, which is installed if missing.

```r
if (!require("devtools", quietly = TRUE))
    install.packages("devtools")
devtools::install_github("SONGDONGYUAN1994/scDesign3")
```

--------------------------------

### Example: perform_lrt with constructed data and fitted marginal models

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/perform_lrt.html

Complete example from the documentation. It constructs two datasets with different correlation structures, fits marginal models with and without a pseudotime smooth term, extracts the fits, and performs the likelihood ratio test. The output includes a progress bar from the second construct_data call.

```r
data(example_sce)
my_data <- construct_data(
sce = example_sce,
assay_use = "counts",
celltype = "cell_type",
pseudotime = "pseudotime",
spatial = NULL,
other_covariates = NULL,
corr_by = "cell_type"
)

my_data2 <- construct_data(
  sce = example_sce,
  assay_use = "counts",
  celltype = "cell_type",
  pseudotime = "pseudotime",
  spatial = NULL,
  other_covariates = NULL,
  corr_by = "pseudotime",
  ncell = 10000
)
#> fitting ...
#> 
  |                                                                          
  |                                                                      |   0%
  |                                                                          
  |======                                                                |   9%
  |                                                                          
  |=============                                                         |  18%
  |                                                                          
  |===================                                                   |  27%
  |                                                                          
  |=========================                                             |  36%
  |                                                                          
  |================================                                      |  45%
  |                                                                          
  |======================================                                |  55%
  |                                                                          
  |=============================================                         |  64%
  |                                                                          
  |===================================================                   |  73%
  |                                                                          
  |=========================================================             |  82%
  |                                                                          
  |================================================================      |  91%
  |                                                                          
  |======================================================================| 100%

my_marginal1 <- fit_marginal(
  data = my_data,
  mu_formula = "1",
  sigma_formula = "1",
  family_use = "nb",
  n_cores = 1,
  usebam = FALSE
)
my_marginal2 <- fit_marginal(
  data = my_data,
  mu_formula = "s(pseudotime, bs = 'cr', k = 10)",
  sigma_formula = "1",
  family_use = "nb",
  n_cores = 1,
  usebam = FALSE
)
my_fit1 <- lapply(my_marginal1, function(x)x$fit)
my_fit2 <- lapply(my_marginal2, function(x)x$fit)
my_pvalue <- perform_lrt(my_fit2, my_fit1)
```

--------------------------------

### Read and inspect example_sce from figshare

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-conditionEffect-vignette.html

Downloads the example PBMC dataset from figshare and prints its structure. The dataset is a SingleCellExperiment with 1000 genes and 13999 cells, containing condition information in colData.

```r
example_sce <- readRDS((url("https://figshare.com/ndownloader/files/40581977")))
print(example_sce)
#> class: SingleCellExperiment 
#> dim: 1000 13999 
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(1000): FTL HBB ... AC022182.3 SCIN
#> rowData names(0):
#> colnames(13999): AAACATACATTTCC.1 AAACATACCAGAAA.1 ... TTTGCATGCTAAGC.1
#>   TTTGCATGGGACGA.1
#> colData names(8): orig.ident nCount_RNA ... cell_type condition
#> reducedDimNames(0):
#> mainExpName: RNA
#> altExpNames(0):
```

--------------------------------

### Fit vine copula with example data

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/fit_copula.html

Example workflow that constructs data, fits marginal models, and then fits a vine copula. The output shows the expected console messages, including the time difference for the vine copula estimation.

```r
data(example_sce)
  my_data <- construct_data(
  sce = example_sce,
  assay_use = "counts",
  celltype = "cell_type",
  pseudotime = "pseudotime",
  spatial = NULL,
  other_covariates = NULL,
  corr_by = "1"
  )
  my_marginal <- fit_marginal(
  data = my_data,
  mu_formula = "s(pseudotime, bs = 'cr', k = 10)",
  sigma_formula = "1",
  family_use = "nb",
  n_cores = 1,
  usebam = FALSE
  )
  my_copula <- fit_copula(
  sce = example_sce,
  assay_use = "counts",
  marginal_list = my_marginal,
  family_use = c(rep("nb", 5), rep("zip", 5)),
  copula = "vine",
  n_cores = 1,
  input_data = my_data$dat
  )
#> Convert Residuals to Uniform
#> Converting End
#> Copula group 1 starts
#> Vine Copula Estimation Starts
#> Time difference of 0.1842949 secs
#> Vine Copula Estimation Ends
```

--------------------------------

### extract_para example with construct_data, fit_marginal, and fit_copula

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/extract_para.html

Example workflow for extract_para. It constructs data, fits marginal models, fits a vine copula, and then calls extract_para with the fitted marginal list and new covariates. The output shows the copula fitting progress messages.

```r
  data(example_sce)
  my_data <- construct_data(
  sce = example_sce,
  assay_use = "counts",
  celltype = "cell_type",
  pseudotime = "pseudotime",
  spatial = NULL,
  other_covariates = NULL,
  corr_by = "1"
  )
  my_marginal <- fit_marginal(
  data = my_data,
  mu_formula = "s(pseudotime, bs = 'cr', k = 10)",
  sigma_formula = "1",
  family_use = "nb",
  n_cores = 1,
  usebam = FALSE
  )
  my_copula <- fit_copula(
  sce = example_sce,
  assay_use = "counts",
  marginal_list = my_marginal,
  family_use = c(rep("nb", 5), rep("zip", 5)),
  copula = "vine",
  n_cores = 1,
  input_data = my_data$dat
  )
#> Convert Residuals to Uniform
#> Converting End
#> Copula group 1 starts
#> Vine Copula Estimation Starts
#> Time difference of 0.1801219 secs
#> Vine Copula Estimation Ends
  my_para <- extract_para(
    sce = example_sce,
    marginal_list = my_marginal,
    n_cores = 1,
    family_use = c(rep("nb", 5), rep("zip", 5)),
    new_covariate = my_data$new_covariate,
    data = my_data$dat
  )
```

--------------------------------

### Example: fitting marginal models with fit_marginal

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/fit_marginal.html

Example usage of fit_marginal with a constructed data object. It uses the 'nb' family, a smooth term for pseudotime, and one core. The example assumes example_sce is available.

```r
data(example_sce)
  my_data <- construct_data(
  sce = example_sce,
  assay_use = "counts",
  celltype = "cell_type",
  pseudotime = "pseudotime",
  spatial = NULL,
  other_covariates = NULL,
  corr_by = "1"
  )
  my_marginal <- fit_marginal(
  data = my_data,
  mu_formula = "s(pseudotime, bs = 'cr', k = 10)",
  sigma_formula = "1",
  family_use = "nb",
  n_cores = 1,
  usebam = FALSE
  )
```

--------------------------------

### Simulate with vine copula using scdesign3

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-copulaCompare-vignette.html

Simulates a new dataset using the vine copula. Uses the same parameters as the Gaussian example but with copula = "vine". The result is stored in Zhengmix4eq_simu_sce_vine.

```R
set.seed(123)
Zhengmix4eq_simu_sce_vine <- scdesign3(sce = Zhengmix4eq_sce_sub,
                            celltype = 'cell_type',
                            pseudotime = NULL,
                            spatial = NULL,
                            other_covariates = NULL,
                            corr_formula = "cell_type",
                            mu_formula = "cell_type",
                            sigma_formula = "cell_type",
                            n_cores = 2,
                            copula = "vine",
                            assay_use = "normcounts",
                            family_use = "nb",
                            pseudo_obs = TRUE, return_model = TRUE)
```

--------------------------------

### Inspect colData with batch information

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-batchEffect-vignette.html

Displays the first six rows of the colData from the example SingleCellExperiment object, showing the batch column and other metadata such as cell type and experiment.

```r
head(colData(example_sce))
#> DataFrame with 6 rows and 13 columns
#>                                        phenoid orig.ident nCount_RNA
#>                                    <character>   <factor>  <numeric>
#> pbmc2_10X_V2_AAACCTGAGATGGGTC           B cell      pbmc2       2360
#> pbmc2_10X_V2_AAACCTGAGCGTAATA           B cell      pbmc2       1888
#> pbmc2_10X_V2_AAACCTGAGCTAGGCA Cytotoxic T cell      pbmc2       3456
#> pbmc2_10X_V2_AAACCTGAGGGTCTCC   Dendritic cell      pbmc2       3802
#> pbmc2_10X_V2_AAACCTGGTCCGAACC      CD4+ T cell      pbmc2       3826
#> pbmc2_10X_V2_AAACCTGTCGTCCGTT      CD4+ T cell      pbmc2       2345
#>                               nFeature_RNA       nGene        nUMI
#>                                  <integer> <character> <character>
#> pbmc2_10X_V2_AAACCTGAGATGGGTC         1044        1044        2360
#> pbmc2_10X_V2_AAACCTGAGCGTAATA          803         803        1888
#> pbmc2_10X_V2_AAACCTGAGCTAGGCA         1372        1372        3456
#> pbmc2_10X_V2_AAACCTGAGGGTCTCC         1519        1519        3802
#> pbmc2_10X_V2_AAACCTGGTCCGAACC         1451        1451        3826
#> pbmc2_10X_V2_AAACCTGTCGTCCGTT          931         931        2345
#>                                     percent.mito     Cluster  Experiment
#>                                      <character> <character> <character>
#> pbmc2_10X_V2_AAACCTGAGATGGGTC 0.0419491525423729           2       pbmc2
#> pbmc2_10X_V2_AAACCTGAGCGTAATA 0.0413135593220339           2       pbmc2
#> pbmc2_10X_V2_AAACCTGAGCTAGGCA 0.0353009259259259           1       pbmc2
#> pbmc2_10X_V2_AAACCTGAGGGTCTCC 0.0420831141504471           6       pbmc2
#> pbmc2_10X_V2_AAACCTGGTCCGAACC 0.0371144798745426           0       pbmc2
#> pbmc2_10X_V2_AAACCTGTCGTCCGTT 0.0652452025586354           0       pbmc2
#>                                          Method    ident             batch
#>                                     <character> <factor>       <character>
#> pbmc2_10X_V2_AAACCTGAGATGGGTC 10x Chromium (v2)    pbmc2 10x Chromium (v2)
#> pbmc2_10X_V2_AAACCTGAGCGTAATA 10x Chromium (v2)    pbmc2 10x Chromium (v2)
#> pbmc2_10X_V2_AAACCTGAGCTAGGCA 10x Chromium (v2)    pbmc2 10x Chromium (v2)
#> pbmc2_10X_V2_AAACCTGAGGGTCTCC 10x Chromium (v2)    pbmc2 10x Chromium (v2)
#> pbmc2_10X_V2_AAACCTGGTCCGAACC 10x Chromium (v2)    pbmc2 10x Chromium (v2)
#> pbmc2_10X_V2_AAACCTGTCGTCCGTT 10x Chromium (v2)    pbmc2 10x Chromium (v2)
#>                                      cell_type
#>                                    <character>
#> pbmc2_10X_V2_AAACCTGAGATGGGTC           B cell
#> pbmc2_10X_V2_AAACCTGAGCGTAATA           B cell
#> pbmc2_10X_V2_AAACCTGAGCTAGGCA Cytotoxic T cell
#> pbmc2_10X_V2_AAACCTGAGGGTCTCC   Dendritic cell
#> pbmc2_10X_V2_AAACCTGGTCCGAACC      CD4+ T cell
#> pbmc2_10X_V2_AAACCTGTCGTCCGTT      CD4+ T cell
```

--------------------------------

### Load example spatial SingleCellExperiment from figshare

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-DEanalysis-vignette.html

Reads a remote RDS file from figshare containing a SingleCellExperiment with 1000 genes and 2696 cells, with counts and logcounts assays and cell_type column. The example output shows the class, dimensions, and assay information.

```R
example_sce <- readRDS((url("https://figshare.com/ndownloader/files/40582019")))
print(example_sce)
#> class: SingleCellExperiment 
#> dim: 1000 2696 
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(1000): Calb2 Gng4 ... Fndc5 Gda
#> rowData names(0):
#> colnames(2696): AAACAAGTATCTCCCA-1 AAACACCAATAACTGC-1 ...
#>   TTGTTTCATCCAGG-1 TTGTTTCCATACAACT-1
#> colData names(12): orig.ident nCount_Spatial ... spatial2 cell_type
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
```

--------------------------------

### Read example SingleCellExperiment data from figshare

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-parallelization-vignette.html

Loads the preprocessed pancreatic endocrinogenesis dataset (top 1000 highly variable genes, single trajectory) from a figshare URL using readRDS and prints it. The output shows the SingleCellExperiment structure with 1000 genes and 2087 cells.

```r
example_sce <- readRDS((url("https://figshare.com/ndownloader/files/40581992")))
print(example_sce)
```

--------------------------------

### Print cell type table from example_sce

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-conditionEffect-vignette.html

Prints the cell type counts in the example dataset. The output shows the number of cells per cell type, with CD14 Mono being the most abundant and Eryth the least.

```r
print(table(colData(example_sce)$cell_type))
#> 
#>    CD14 Mono  CD4 Naive T CD4 Memory T    CD16 Mono            B        CD8 T 
#>         4362         2504         1762         1044          978          814 
#>  T activated           NK           DC  B Activated           Mk          pDC 
#>          633          619          472          388          236          132 
#>        Eryth 
#>           55
```

--------------------------------

### Construct data with construct_data()

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-DEanalysis-vignette.html

Initializes the data structure for simulation using the example_sce object. Sets a random seed for reproducibility.

```R
set.seed(123)
example_data <- construct_data(
    sce = example_sce,
    assay_use = "counts",
    celltype = "cell_type",
    pseudotime = NULL,
    spatial = NULL,
    other_covariates = NULL,
    corr_by = "1"
  )
```

--------------------------------

### Load required libraries

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-spatial-deconvolution.html

Loads scDesign3, SingleCellExperiment, ggplot2, dplyr, viridis, IOBR, and scatterpie for the tutorial.

```R
library(scDesign3)
library(SingleCellExperiment)
library(ggplot2)
library(dplyr)
library(viridis)
library(IOBR)
library(scatterpie)
theme_set(theme_bw())
```

--------------------------------

### Load required libraries for scDesign3 DE benchmarking

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-DEanalysis-vignette.html

Loads all R packages needed for the tutorial, including scDesign3, SingleCellExperiment, Seurat, scran, DESeq2, tradeSeq, and spatial DE tools. Run this before any other code in the vignette.

```R
library(scDesign3)
library(SingleCellExperiment)
library(ggplot2)
library(DuoClustering2018)
library(Seurat)
library(SeuratObject)
library(scran)
library(parallel)
library(monocle3)
library(DESeq2)
library(BiocParallel)
library(NBAMSeq)
library(PseudotimeDE)
library(dplyr)
library(tradeSeq)
library(reshape2)
library(spatialDE)
library(SPARK)
library(scales)
library(tidyr)
theme_set(theme_bw())
```

--------------------------------

### Subset to top 100 genes

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-parallelization-vignette.html

Reduces the example SingleCellExperiment to the first 100 genes to save computational time before running scdesign3.

```R
example_sce <- example_sce[1:100, ]
```

--------------------------------

### Subset to B cells and regulatory T cells

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-DEanalysis-vignette.html

Filters the example SingleCellExperiment to only B cells and regulatory T cells, and adds a 'cell_type' factor column to the colData for downstream analysis.

```R
selected_cells <- which(colData(example_sce)$phenoid %in% c("b.cells","regulatory.t"))
example_sce <- example_sce[,selected_cells]
colData(example_sce)$cell_type <- as.factor(colData(example_sce)$phenoid)
example_sce
head(colData(example_sce))
```

--------------------------------

### Run scdesign3 with default parallelization (mcmapply)

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-parallelization-vignette.html

Uses the default parallelization from the 'parallel' package. No need to set the parallelization parameter.

```R
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",
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    DT = TRUE,
    pseudo_obs = FALSE,
    return_model = FALSE,
    nonzerovar = FALSE
  )
```

--------------------------------

### Run scdesign3 with example_sce and Gaussian copula

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/scdesign3.html

Runs scdesign3 on the example_sce dataset with a negative binomial family, a pseudotime smooth term for mu, and a Gaussian copula. The output shows the progress messages and the sequence of fitting steps, from input data construction through marginal and copula fitting to new data generation.

```r
data(example_sce)
my_simu <- scdesign3(
sce = example_sce,
assay_use = "counts",
celltype = "cell_type",
pseudotime = "pseudotime",
spatial = NULL,
other_covariates = NULL,
mu_formula = "s(pseudotime, bs = 'cr', k = 10)",
sigma_formula = "1",
family_use = "nb",
n_cores = 2,
usebam = FALSE,
edf_flexible = FALSE,
corr_formula = "pseudotime",
copula = "gaussian",
if_sparse = TRUE,
DT = TRUE,
pseudo_obs = FALSE,
ncell = 1000,
return_model = FALSE
)
#> Input Data Construction Start
#> fitting ...
#> 
  |                                                                          
  |                                                                      |   0%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |======================================================================| 100%
#> Input Data Construction End
#> Start Marginal Fitting
#> Marginal Fitting End
#> Start Copula Fitting
#> Convert Residuals to Multivariate Gaussian
#> Converting End
#> Copula group 1 starts
#> Copula group 2 starts
#> Copula group 3 starts
#> Copula group 4 starts
#> Copula Fitting End
#> Start Parameter Extraction
#> Parameter
#> Extraction End
#> Start Generate New Data
#> Use Copula to sample a multivariate quantile matrix
#> Sample Copula group 1 starts
#> Sample Copula group 2 starts
#> Sample Copula group 3 starts
#> Sample Copula group 4 starts
#> New Data Generating End
```

--------------------------------

### Construct input dataset for step-by-step simulation

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-cellType-vignette.html

Uses construct_data to prepare the input dataset for step-by-step simulation, specifying assay_use='counts', celltype='cell_type', and corr_by='cell_type'.

```r
example_data <- construct_data(
    sce = sce,
    assay_use = "counts",
    celltype = "cell_type",
    pseudotime = NULL,
    spatial = NULL,
    other_covariates = NULL,
    corr_by = "cell_type"
  )
```

--------------------------------

### Read and print reference scATAC-seq data

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-scATACseq-vignette.html

Downloads the example SingleCellExperiment object from figshare and prints its structure. The object contains 1133 peaks and 7034 cells with assays 'counts' and 'logcounts'.

```R
example_sce <- readRDS((url("https://figshare.com/ndownloader/files/40581962")))
print(example_sce)
#> class: SingleCellExperiment 
#> dim: 1133 7034 
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(1133): chr6-44025105-44028184 chr2-113581628-113594911 ...
#>   chr2-233174300-233175622 chr3-13015237-13015864
#> rowData names(0):
#> colnames(7034): AAACGAAAGAGCGAAA-1 AAACGAAAGAGTTTGA-1 ...
#>   TTTGTGTTCTACTTTG-1 TTTGTGTTCTTGTGCC-1
#> colData names(31): orig.ident nCount_peaks ... ident cell_type
#> reducedDimNames(2): LSI UMAP
#> mainExpName: peaks
#> altExpNames(0)
```

--------------------------------

### Inspect colData of example_sce

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-conditionEffect-vignette.html

Displays the first six rows of the colData, showing the condition column with values CTRL and the cell type annotations for the subsetted dataset.

```r
head(colData(example_sce))
#> DataFrame with 6 rows and 8 columns
#>                   orig.ident nCount_RNA nFeature_RNA        stim
#>                  <character>  <numeric>    <integer> <character>
#> AAACATACATTTCC.1 IMMUNE_CTRL       3017          248        CTRL
#> AAACATACCAGAAA.1 IMMUNE_CTRL       2481          713        CTRL
#> AAACATACCTCGCT.1 IMMUNE_CTRL       3420          850        CTRL
#> AAACATACGGCATT.1 IMMUNE_CTRL       1581          557        CTRL
#> AAACATTGCTTCGC.1 IMMUNE_CTRL       2536          669        CTRL
#> AAACGCACTCGCCT.1 IMMUNE_CTRL       3563          908        CTRL
#>                  seurat_annotations       ident cell_type   condition
#>                            <factor>    <factor>  <factor> <character>
#> AAACATACATTTCC.1          CD14 Mono IMMUNE_CTRL CD14 Mono        CTRL
#> AAACATACCAGAAA.1          CD14 Mono IMMUNE_CTRL CD14 Mono        CTRL
#> AAACATACCTCGCT.1          CD14 Mono IMMUNE_CTRL CD14 Mono        CTRL
#> AAACATACGGCATT.1          CD14 Mono IMMUNE_CTRL CD14 Mono        CTRL
#> AAACATTGCTTCGC.1          CD14 Mono IMMUNE_CTRL CD14 Mono        CTRL
#> AAACGCACTCGCCT.1          CD14 Mono IMMUNE_CTRL CD14 Mono        CTRL
```

--------------------------------

### Run scdesign3 with pbmcapply parallelization

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-parallelization-vignette.html

Sets parallelization to 'pbmcmapply' to show a progress bar during parallel execution.

```R
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",
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    DT = TRUE,
    pseudo_obs = FALSE,
    return_model = FALSE,
    nonzerovar = FALSE,
    parallelization = "pbmcmapply"
  )
```

--------------------------------

### Subset example_sce to top genes and two cell types

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-conditionEffect-vignette.html

Subsets the example dataset to the top 100 genes and only the CD14 Mono and B cell types, then removes unused cell type factor levels to save computational time.

```r
example_sce <- example_sce[1:100, colData(example_sce)$cell_type %in% c("CD14 Mono", "B")]

## Remove unused cell type levels
colData(example_sce)$cell_type <- droplevels(colData(example_sce)$cell_type)
```

--------------------------------

### Load required R packages for pseudotime GOF evaluation

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-pseudotimeGOF-vignette.html

Loads the required R packages for the tutorial: scDesign3, dyngen, SingleCellExperiment, ggplot2, and dplyr. Sets the ggplot2 theme to black and white.

```r
library(scDesign3)
library(dyngen)
library(SingleCellExperiment)
library(ggplot2)
library(dplyr)
theme_set(theme_bw())
```

--------------------------------

### Subset Zhengmix4eq to top HVGs and TFs

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-copulaCompare-vignette.html

Selects the top 100 highly variable genes, adds human transcription factors, and filters out housekeeping genes and genes starting with RP or TMSB. Subsets the SingleCellExperiment to these genes and prints the resulting object, showing 119 genes and 3555 cells.

```r
humantfs <- read_csv("http://humantfs.ccbr.utoronto.ca/download/v_1.01/TF_names_v_1.01.txt", col_names = FALSE)
stats <- modelGeneVar(Zhengmix4eq_sce)
chosen <- getTopHVGs(stats, n = 100)
### Filter out some HKGs.
chosen <- union(chosen, humantfs$X1)
chosen <- chosen[!stringr::str_starts(chosen, "RP")]
chosen <- chosen[!stringr::str_starts(chosen, "TMSB")]
chosen <- chosen[!chosen %in% c("B2M", "MALAT1", "ACTB", "ACTG1", "GAPDH", "FTL", "FTH1")]
Zhengmix4eq_sce_sub <- Zhengmix4eq_sce[rownames(Zhengmix4eq_sce) %in% chosen, ]
print(Zhengmix4eq_sce_sub)
#> class: SingleCellExperiment 
#> dim: 119 3555 
#> metadata(1): log.exprs.offset
#> assays(3): counts logcounts normcounts
#> rownames(119): TPT1 CD74 ... PHF1 IRF7
#> rowData names(10): id ... total_counts log10_total_counts
#> colnames(3555): b.cells6276 b.cells6144 ... regulatory.t1084
#>   regulatory.t9696
#> colData names(16): dataset barcode ... sizeFactor cell_type
#> reducedDimNames(2): PCA TSNE
#> mainExpName: NULL
#> altExpNames(0):
```

--------------------------------

### Read and inspect reference spatial data (example_sce)

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-spatial-vignette.html

Downloads and reads a SingleCellExperiment object from figshare containing Visium spatial transcriptomics data with 1000 genes and 2696 cells. Prints the object to show its structure.

```R
example_sce <- readRDS((url("https://figshare.com/ndownloader/files/40582019")))
print(example_sce)
```

--------------------------------

### Initialize q-values matrix for DE methods

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-DEanalysis-vignette.html

Creates a matrix to store BH-corrected p-values from four DE methods: Monocle3-DE, tradeSeq, NBAMseq, and PseudotimeDE.

```r
qvals <- matrix(0, ncol = 4, nrow = dim(simu_sce)[1])
colnames(qvals) <- c("Monocle3-DE","tradeSeq","NBAMseq","PseudotimeDE")
rownames(qvals) <- rownames(simu_sce)
```

--------------------------------

### Construct synthetic spatial data with construct_data()

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-DEanalysis-vignette.html

Builds a data structure for spatial simulation using the cell type and spatial coordinates from example_sce. Sets the seed to 1 for reproducibility and uses "counts" as the assay.

```R
set.seed(1)
example_data <- construct_data(
    sce = example_sce,
    assay_use = "counts",
    celltype = "cell_type",
    pseudotime = NULL,
    spatial = c("spatial1", "spatial2"),
    other_covariates = NULL,
    corr_by = "1"
  )
```

--------------------------------

### Run scdesign3 with BiocParallel using SnowParam

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-parallelization-vignette.html

Uses pbmapply from BiocParallel with SnowParam for parallelization.

```R
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",
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    DT = TRUE,
    pseudo_obs = FALSE,
    return_model = FALSE,
    nonzerovar = FALSE,
    parallelization = "pbmapply",
    BPPARAM = SnowParam()
  )
```

--------------------------------

### Load required R packages for copula comparison

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-copulaCompare-vignette.html

Loads the required R packages for the tutorial, including scDesign3, DuoClustering2018, scater, tidygraph, ggrepel, scran, dplyr, igraph, tidyverse, and rvinecopulib. Sets the default ggplot2 theme to black and white.

```r
library(scDesign3)
library(DuoClustering2018)
library(scater)
library(tidygraph)
library(ggrepel)
library(scran)
library(dplyr)
library(igraph)
library(tidyverse)
library(rvinecopulib)
theme_set(theme_bw())
```

--------------------------------

### Generate reference dataset with dyngen and store ground truth pseudotime

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-pseudotimeGOF-vignette.html

Generates a reference dataset with ground truth pseudotime using dyngen. Sets a random seed for reproducibility, creates a linear backbone, initializes the model with 500 cells and specified numbers of transcription factors, target genes, and housekeeping genes, then generates a SingleCellExperiment object and stores the true pseudotime in its colData.

```r
set.seed(123)
backbone <- backbone_linear_simple()
config <-
  initialise(
    backbone = backbone,
    num_cells = 500,
    num_tfs = nrow(backbone$module_info),
    num_targets = 100,
    num_hks = 50,
    verbose = FALSE
  )
out <- generate_dataset(
  config,
  format = "sce",
  make_plots = FALSE
  )
example_sce <- out$dataset
colData(example_sce)$pseudotime <- out$model$experiment$cell_info$time
```

--------------------------------

### Load scDesign3 and ggplot2 libraries

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-introduction-vignette.html

Loads the required packages and sets a black-and-white theme for plots.

```R
library(scDesign3)
library(ggplot2)
theme_set(theme_bw())
```

--------------------------------

### Run scdesign3 with BiocParallel using MulticoreParam

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-parallelization-vignette.html

Uses pbmapply from BiocParallel with MulticoreParam for parallelization.

```R
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",
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    DT = TRUE,
    pseudo_obs = FALSE,
    return_model = FALSE,
    nonzerovar = FALSE,
    parallelization = "pbmapply",
    BPPARAM = MulticoreParam()
  )
```

--------------------------------

### Load example_sce dataset

Source: https://songdongyuan1994.github.io/scDesign3/docs/reference/example_sce.html

Loads the example_sce dataset, a SingleCellExperiment object containing both cell type and pseudotime information. The dataset has 10 rows (genes) and 1289 columns (cells).

```r
data("example_sce")
```

--------------------------------

### Show session information

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-librarySize-vignette.html

Display the R session information, including version, platform, and locale.

```R
sessionInfo()
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.6 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.9.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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/Los_Angeles
#> tzcode source: system (glibc)
#> 
```

--------------------------------

### Run scDesign3 simulation with spatial coordinates

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-spatial-vignette.html

Calls scdesign3() to simulate new spatial data. Uses a Gaussian process smooth term for spatial coordinates in mu_formula, negative binomial family, and parallelization with pbmcapply. Sets a seed for reproducibility.

```R
set.seed(123)
example_simu <- scdesign3(
    sce = example_sce,
    assay_use = "counts",
    celltype = "cell_type",
    pseudotime = NULL,
    spatial = c("spatial1", "spatial2"),
    other_covariates = NULL,
    mu_formula = "s(spatial1, spatial2, bs = 'gp', k= 400)",
    sigma_formula = "1",
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    DT = TRUE,
    pseudo_obs = FALSE,
    return_model = FALSE,
    nonzerovar = FALSE,
    parallelization = "pbmcapply"
  )
```

--------------------------------

### Load required packages and set theme

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-librarySize-vignette.html

Load required R packages and set the default ggplot2 theme to black and white.

```R
library(scDesign3)
library(SingleCellExperiment)
library(ggplot2)
library(DuoClustering2018)
library(scran)
library(tidyverse)
theme_set(theme_bw())
```

--------------------------------

### Load required libraries for scDesign3 spatial simulation

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-spatial-vignette.html

Loads the scDesign3, SingleCellExperiment, ggplot2, dplyr, and viridis packages and sets the ggplot2 theme to theme_bw(). Required for the tutorial.

```R
library(scDesign3)
library(SingleCellExperiment)
library(ggplot2)
library(dplyr)
library(viridis)
theme_set(theme_bw())
```

--------------------------------

### Prepare training data with top 100 HVGs and PCAKmeans results

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-clusterGOF-vignette.html

Extracts the PCAKmeans clustering results, selects the top 100 highly variable genes, and splits the data into a training set of half the cells with a fixed seed.

```r
kmeans_res <- Filter(function(x){all(x$method == "PCAKmeans")}, res_sub_list)
ncell <- ncol(Zhengmix4eq_sce)
ngene <- 100
zheng_sce <- modelGeneVar(Zhengmix4eq_sce)
chosen <- getTopHVGs(zheng_sce, n = ngene)
sce <- Zhengmix4eq_sce[chosen, ]
ntrain <- round(ncell/2)
  
set.seed(123)
train_index <- sample(seq_len(ncell), ntrain, replace = FALSE)
train_sce <- sce[, train_index]
```

--------------------------------

### Simulate data with scDesign3 using library size offset

Source: https://songdongyuan1994.github.io/scDesign3/docs/articles/scDesign3-librarySize-vignette.html

Run scDesign3 simulation with library size as an offset in the mu formula. The library size is log-transformed because the link function for mu of the negative binomial distribution in GAMLSS is log.

```R
set.seed(123)
example_simu <- scdesign3(
    sce = sce,
    assay_use = "counts",
    celltype = "cell_type",
    pseudotime = NULL,
    spatial = NULL,
    other_covariates = "library",
    mu_formula = "cell_type + offset(log(library))",
    sigma_formula = "1",
    family_use = "nb",
    n_cores = 2,
    usebam = FALSE,
    corr_formula = "1",
    copula = "gaussian",
    DT = TRUE,
    pseudo_obs = FALSE,
    return_model = FALSE,
    nonzerovar = FALSE,
    parallelization = "pbmcmapply",
    important_feature = "auto"
  )
```
Tokens:
Raw