Inputs and Outputs

library(Seurat)
#> Loading required package: methods
#> Loading required package: SeuratObject
#> Loading required package: sp
#> 
#> Attaching package: 'SeuratObject'
#> The following object is masked from 'package:base':
#> 
#>     intersect
library(curl)
#> Using libcurl 7.81.0 with OpenSSL/3.0.2
library(SimBu)

This chapter covers the different input and output options of the package in detail.

Input

The input for your simulations is always a SummarizedExperiment object. You can create this object with different constructing functions, which will explained below. It is also possible to merge multiple datasets objects into one.
Sfaira is not covered in this vignette, but in “Public Data Integration”.

Custom data

Using existing count matrices and annotations is already covered in the “Getting started” vignette; this section will explain some minor details.

When generating a dataset with your own data, you need to provide the count_matrix parameter of dataset(); additionally you can provide a TPM matrix with the tpm_matrix. This will then lead to two simulations, one based on counts and one based on TPMs. For either of them, genes are located in the rows, cells in the columns.
Additionally, an annotation table is needed, with the cell-type annotations. It needs to consist of at least out of 2 columns: ID and cell_type, where ID has to be identical to the column names of the provides matrix/matrices. If not all cells appear in the annotation or matrix, the intersection of both is used to generate the dataset.

Here is some example data:

counts <- Matrix::Matrix(matrix(stats::rpois(3e5, 5), ncol = 300), sparse = TRUE)
tpm <- Matrix::Matrix(matrix(stats::rpois(3e5, 5), ncol = 300), sparse = TRUE)
tpm <- Matrix::t(1e6 * Matrix::t(tpm) / Matrix::colSums(tpm))
colnames(counts) <- paste0("cell-", rep(1:300))
colnames(tpm) <- paste0("cell-", rep(1:300))
rownames(counts) <- paste0("gene-", rep(1:1000))
rownames(tpm) <- paste0("gene-", rep(1:1000))
annotation <- data.frame(
  "ID" = paste0("cell-", rep(1:300)),
  "cell_type" = c(
    rep("T cells CD4", 50),
    rep("T cells CD8", 50),
    rep("Macrophages", 100),
    rep("NK cells", 10),
    rep("B cells", 70),
    rep("Monocytes", 20)
  ),
  row.names = paste0("cell-", rep(1:300))
)
seurat_obj <- Seurat::CreateSeuratObject(counts = counts, assay = "gene_expression", meta.data = annotation)
# store normalized matrix in the 'data' layer
SeuratObject::LayerData(seurat_obj, assay = "gene_expression", layer = "data") <- tpm
seurat_obj
#> An object of class Seurat 
#> 1000 features across 300 samples within 1 assay 
#> Active assay: gene_expression (1000 features, 0 variable features)
#>  2 layers present: counts, data

Seurat

It is possible to use a Seurat object to build a dataset; give the name of the assay containing count data in the counts slot, the name of the column in the meta table with the unique cell IDs and the name of the column in the meta table with the cell type identifier. Additionally you may give the name of the assay containing TPM data in the counts slot.

ds_seurat <- SimBu::dataset_seurat(
  seurat_obj = seurat_obj,
  counts_layer = "counts",
  cell_id_col = "ID",
  cell_type_col = "cell_type",
  tpm_layer = "data",
  name = "seurat_dataset"
)
#> Filtering genes...
#> Created dataset.

h5ad files

It is possible to use an h5ad file directly, a file format which stores AnnData objects. As h5ad files can store cell specific information in the obs layer, no additional annotation input for SimBu is needed.
Note: if you want both counts and tpm data as input, you will have to provide two files; the cell annotation has to match between these two files. As SimBu expects the cells to be in the columns and genes/features in the rows of the input matrix, but this is not necessarily the case for anndata objects https://falexwolf.de/img/scanpy/anndata.svg, SimBu can handle h5ad files with cells in the obs or var layer. If your cells are in obs, use cells_in_obs=TRUE and FALSE otherwise. This will also automatically transpose the matrix. To know, which columns in the cell annotation layer correspond to the cell identifiers and cell type labels, use the cell_id_col and cell_type_col parameters, respectively.

As this function uses the SimBu python environment to read the h5ad files and extract the data, it may take some more time to initialize the conda environment at the first usage only.

# example h5ad file, where cell type info is stored in `obs` layer
# h5 <- system.file("extdata", "anndata.h5ad", package = "SimBu")
# ds_h5ad <- SimBu::dataset_h5ad(
#  h5ad_file_counts = h5,
#  name = "h5ad_dataset",
#  cell_id_col = 0, # this will use the rownames of the metadata as cell identifiers
#  cell_type_col = "group", # this will use the 'group' column of the metadata as cell type info
#  cells_in_obs = TRUE # in case your cell information is stored in the var layer, switch to FALSE
# )

Merging datasets

You are able to merge multiple datasets by using the dataset_merge function:

ds <- SimBu::dataset(
  annotation = annotation,
  count_matrix = counts,
  tpm_matrix = tpm,
  name = "test_dataset"
)
#> Filtering genes...
#> Created dataset.
ds_multiple <- SimBu::dataset_merge(
  dataset_list = list(ds_seurat, ds),
  name = "ds_multiple"
)
#> Filtering genes...
#> Created dataset.

Output

The simulation object contains three named entries:

bulk: a SummarizedExperiment object with the pseudo-bulk dataset(s) stored in the assays. They can be accessed like this:

simulation <- SimBu::simulate_bulk(
  data = ds_multiple,
  scenario = "random",
  scaling_factor = "NONE",
  nsamples = 10, ncells = 100
)
#> Finished simulation.
dim(SummarizedExperiment::assays(simulation$bulk)[["bulk_counts"]])
#> [1] 1000   10
dim(SummarizedExperiment::assays(simulation$bulk)[["bulk_tpm"]])
#> [1] 1000   10

If only the count matrix was given to the dataset initially, only the bulk_counts assay is filled.

cell_fractions: a table where rows represent the simulated samples and columns represent the different simulated cell-types. The entries in the table store the specific cell-type fraction per sample.
scaling_vector: a named list, with the used scaling value for each cell from the single cell dataset.