ESCO

Here ESCO method will be demonstrated clearly and hope that this document can help you.

Estimating parameters from a real dataset

Before simulating datasets, it is important to estimate some essential parameters from a real dataset in order to make the simulated data more real.

library(simmethods)
library(SingleCellExperiment)
# Load data
ref_data <- simmethods::data
dim(ref_data)
# [1] 4000  160

Using simmethods::ESCO_estimation command to execute the estimation step.

estimate_result <- simmethods::ESCO_estimation(
  ref_data = ref_data,
  verbose = T,
  seed = 10
)
# Registered S3 method overwritten by 'DescTools':
#   method         from  
#   reorder.factor gplots
# Registered S3 methods overwritten by 'registry':
#   method               from 
#   print.registry_field proxy
#   print.registry_entry proxy
# Estimating parameters using ESCO

ESCO is not stable, and some datasets can not be estimated due to the failing estimation.

Simulating datasets using ESCO

After estimating parameter from a real dataset, we will simulate a dataset based on the learned parameters with different scenarios.

  1. Datasets with default parameters
  2. Determin the number of cells and genes
  3. Simulate two or more groups

Datasets with default parameters

The reference data contains 160 cells and 4000 genes, if we simulate datasets with default parameters and then we will obtain a new data which has the same size as the reference data. In addtion, the simulated dataset will have one group of cells.

simulate_result <- simmethods::ESCO_simulation(
  parameters = estimate_result[["estimate_result"]],
  return_format = "SCE",
  seed = 111
)
# nCells: 160
# nGenes: 4000
# nGroups: 1
# de.group: 0.1
SCE_result <- simulate_result[["simulate_result"]]
dim(SCE_result)
# [1] 4000  160
head(colData(SCE_result))
# DataFrame with 6 rows and 2 columns
#         cell_name       group
#       <character> <character>
# Cell1       Cell1      Group1
# Cell2       Cell2      Group1
# Cell3       Cell3      Group1
# Cell4       Cell4      Group1
# Cell5       Cell5      Group1
# Cell6       Cell6      Group1

Determin the number of cells and genes

In ESCO, we can set nCells and nGenes to specify the number of cells and genes.

Here, we simulate a new dataset with 500 cells and 1000 genes:

simulate_result <- simmethods::ESCO_simulation(
  parameters = estimate_result[["estimate_result"]],
  return_format = "list",
  other_prior = list(nCells = 500,
                     nGenes = 1000),
  seed = 111
)
# nCells: 500
# nGenes: 1000
# nGroups: 1
# de.group: 0.1
result <- simulate_result[["simulate_result"]][["count_data"]]
dim(result)
# [1] 1000  500

Simulate two or more groups

In ESCO, we can not set nGroups directly and should set prob.group instead. For example, if we want to simulate 2 groups, we can type other_prior = list(prob.group = c(0.5, 0.5)). Note that the sum of prob.group numeric vector must equal to 1, so we can also set prob.group = c(0.3, 0.7).

In addtion, if we want to simulate three or more groups, we should obey the rules:

  • The length of prob.group vector must always equal to the number of groups.
  • The sum of prob.group numeric vector must equal to 1.

For demonstration, we will simulate three groups using the learned parameters.

simulate_result <- simmethods::ESCO_simulation(
  parameters = estimate_result[["estimate_result"]],
  return_format = "list",
  other_prior = list(nCells = 500,
                     nGenes = 1000,
                     prob.group = c(0.1, 0.3, 0.6)),
  seed = 111
)
# nCells: 500
# nGenes: 1000
# nGroups: 3
# de.group: 0.1
result <- simulate_result[["simulate_result"]][["count_data"]]
dim(result)
# [1] 1000  500
## cell information
cell_info <- simulate_result[["simulate_result"]][["col_meta"]]
table(cell_info$group)
# 
# Group1 Group2 Group3 
#     47    168    285
## gene information
gene_info <- simulate_result[["simulate_result"]][["row_meta"]]
### the proportion of DEGs
table(gene_info$de_gene)[2]/nrow(result) ## de.prob = 0.1
# yes 
# 0.1

We can see that the proportion of differential expressed genes is 0.1 (equals to the default). Next, if we want to know the fold change between two groups, we will do division with the groups that we are interested in.

## fc between group2 and group1
fc_group1_to_group2 <- gene_info$DEFacGroup2/gene_info$DEFacGroup1
## fc between group3 and group1
fc_group1_to_group3 <- gene_info$DEFacGroup3/gene_info$DEFacGroup1
## fc between group3 and group2
fc_group2_to_group3 <- gene_info$DEFacGroup3/gene_info$DEFacGroup2