---
title: "2. Simulating Multimodal Single-cell Datasets"
output:
  BiocStyle::html_document:
  toc: true
  toc_depth: 2
vignette: >
  %\VignetteEngine{knitr::knitr}
  %\VignetteIndexEntry{2. Simulating Multimodal Single-cell Datasets}
  %\usepackage[UTF-8]{inputenc}
---
```{r "setup", include=FALSE}
require("knitr")
opts_chunk$set(fig.width=4, fig.height=3)
```

```{r install-packages, include=FALSE, message=FALSE, warning=FALSE, eval=FALSE}
# The following chunk will install all the required packages.
(function() {
  installed <- installed.packages()[,"Package"]
  install <- function(list, fn) {
    pkg <- setdiff(list, installed)
    if (length(pkg)) fn(pkg, dependencies=TRUE)
  }

  r_packages <- c(
    "devtools", "dplyr", "ggplot2", "Rtsne", "rlang",
    "reshape", "ape", "phytools", "repr", "KernelKnn",
    "gridExtra", "parallel", 'foreach', 'phytools', "doParallel",
    "zeallot", "gtools", "gplots", "roxygen2", "usethis"
  )
  install(r_packages, install.packages)

  if (requireNamespace("BiocManager", quietly = TRUE)) {
    bioc_packages <- c('Biobase')
    install(bioc_packages, BiocManager::install)
  }
})()
```

In this tutorial, we will demonstrate how to use scMultiSim to simulate multi-omics data
with different biological effects, including:

- Simulating true RNA counts and ATAC-seq data
- Controlling the cell population and GRN effects
- Adding technical variation and batch effect to the true counts
- Adjusting the parameters to control different biological effects

We first load the package:

```{r load-package, quietly=TRUE, message=FALSE, warning=FALSE}
library("scMultiSim")
```

# Simulating True Counts

scMultiSim first generates the true RNA counts, and then add technical variation and batch effect to the true counts.
To simulate true counts, call `sim_true_counts(options)` where `options` is a
list. You can use `scmultisim_help()` to get help on the options,
or like `scmulti_help("num.cells")` to get help on the options for a specific function.

```{r scmultisim-help, echo = TRUE, results = "hide"}
scmultisim_help("options")
```

## GRN and Differentiation Tree

Before start, we define a utility function to modify a list.
```{r load-dplyr, quietly=TRUE, message=FALSE, warning=FALSE}
library(dplyr)
```
```{r define-list-modify}
list_modify <- function (curr_list, ...) {
  args <- list(...)
  for (i in names(args)) {
    curr_list[[i]] <- args[[i]]
  }
  curr_list
}
```

The minimal input to scMultiSim is a **differentiation tree**, and you can optionally provide
ground truth for GRN and cell-cell interactions.
The differentiation tree is an R phylo object, which can be created using e.g.
`ape::read.tree()` or `ape::rtree()`.
It controls the cell population structure: each node of the tree should represent a cell type,
and connected nodes indicate the differentiation relationship between cell types.
_scMultiSim provides this explicit control on the cell population structure
while preserving all other effects (such as GRN and Cell-Cell Interactions)_,
so you can generate any cell trajectory or clustering structure you want, which is especially useful
for benchmarking trajectory inference and clustering methods.

If generating a continuous population, this tree
specifies the cell differentiation trajectory; if generating a discrete population, the
tips of this tree will be the clusters (cell types are the terminal cell states).

scMultiSim also provides three differentiation trees.
`Phyla5()` and `Phyla3()` return bifurcating trees with 5 and 3 leaves respectively.
`Phyla1()` returns only a single branch, which can be useful when we don't want any specific trajectory.
```{r plot-tree, fig.width = 8, fig.height = 4}
par(mfrow=c(1,2))
Phyla5(plotting = TRUE)
Phyla3(plotting = TRUE)

# It's not possible to plot Phyla1() because it only contains 1 branch connecting two nodes.
Phyla1()
```

If you only need `n` cell clusters without any specific trajectory, you can use code like below to generate a simple tree with `n` leaves.
```{r random-tree}
# tree with four leaves
ape::read.tree(text = "(A:1,B:1,C:1,D:1);")
```

The GRN should be a data frame with 3 columns, each representing the `target`, `regulator`, and `effect`.
The target and regulator should be gene names, which can be integers or strings.
The effect should be a numeric value, indicating the effect of the regulator on the target.

scMultiSim provides two sample GRNs, `GRN_params_100` and `GRN_params_1139`,
which contain 100 and 1139 genes respectively.
Let's load them first.
```{r load-grn}
data(GRN_params_100)
GRN_params <- GRN_params_100
head(GRN_params)
```

## Simulating True Counts

Now, we create the options list for the simulation session.
In the following example, we simulate 500 cells with 50 CIFs.

The number of genes is determined by the option `num.genes` or the number of genes in the GRN.
If `num.genes` is not specified, the number of genes will be the number of unique genes in the GRN,
plus a fraction of genes that are not regulated by any other genes.
this is controlled by the option `unregulated.gene.ratio` (default is 0.1).
Since our `GRN_params` contains 100 gene names, 10% more genes will be added to the simulation,
and the number of genes in the simulated data will be 110.
If you don't need to simulate GRN effects, simply set `GRN = NA`.

The `cif.sigma` controls the variance of the CIFs. Usually, with `cif.sigma` = 0.1,
the trajectory will be very clear, while with `cif.sigma` = 1, the trajectory will be more
noisy. We use `cif.sigma` = 0.5 in this example.

We also have `do.velocity` option to use the Kinetic model to simulate RNA velocity data.

```{r define-options}
set.seed(42)

options <- list(
  GRN = GRN_params,
  num.cells = 300,
  num.cifs = 20,
  cif.sigma = 1,
  tree = Phyla5(),
  diff.cif.fraction = 0.8,
  do.velocity = TRUE
)
```

### Omitting the GRN

Note that the minimal input to scMultiSim is the cell population structure (differentiation tree) and number of cells.
You can omit the GRN by using `GRN = NA`:
```
options <- list(
  GRN = NA
  num.cells = 1000,
  num.genes = 500,
  tree = Phyla5(),
)
```

### Running the Simulation

Now we run the simulation and check what kind of data is in the returned result:
```{r run-simulation}
results <- sim_true_counts(options)
names(results)
```

## Accessing the Results

The return value will be a `scMultiSim Environment` object,
and you can access various data and parameters using the `$` operator.

- `counts`: Gene-by-cell scRNA-seq counts.
- `atac_counts`: Region-by-cell scATAC-seq counts.
- `region_to_gene`: Region-by-gene 0-1 marix indicating the corresponding relationship between chtomatin regions and genes.
- `atacseq_data`: The "clean" scATAC-seq counts without added intrinsic noise.
- `cell_meta`: A dataframe containing cell type labels and pseudotime information.
- `cif`: The CIF used during the simulation.
- `giv`: The GIV used during the simulation.
- `kinetic_params`: The kinetic parameters used during the simulation.
- `.grn`: The GRN used during the simulation.
- `.grn$regulators`: The list of TFs used by all gene-by-TF matrices.
- `.grn$geff`: Gene-by-TF matrix representing the GRN used during the simulation.
- `.n`: Other metadata, e.g. `.n$cells` is the number of cells.

If `do.velocity` is enabled, it has these additional fields:

- `unspliced_counts`: Gene-by-cell unspliced RNA counts.
- `velocity`: Gene-by-cell RNA velocity ground truth.
- `cell_time`: The pseudotime at which the cell counts were generated.

If dynamic GRN is enabled, it has these additional fields:

- `cell_specific_grn`: A list of length `n_cells`. Each element is a gene-by-TF matrix, indicating the cell's GRN.

If cell-cell interaction is enabled, it has these additional fields:

- `grid`: The grid object used during the simulation.
- `grid$get_neighbours(i)`: Get the neighbour cells of cell `i`.
- `cci_locs`: A dataframe containing the X and Y coordinates of each cell.
- `cci_cell_type_param`: A dataframe containing the CCI network ground truth: all ligand-receptor pairs between each pair of cell types.
- `cci_cell_types`: For continuous cell population, the sub-divided cell types along the trajectory used when simulating CCI.

If it is a debug session (`debug = TRUE`), a `sim` field is available,
which is an environment contains all internal states and data structures.

## Visualizing the Results

We can visualize the true counts and ATAC-seq data using `plot_tsne()`:
```{r plot-counts, fig.width = 4, fig.height = 3.5, out.width = "60%"}
plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'True RNA Counts Tsne')
plot_tsne(log2(results$atacseq_data + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'True ATAC-seq Tsne')
```


Since we also have RNA velocity enabled, the `results` also contains the following data:
- `velocity`: the true RNA velocity (genes x cells)
- `unspliced_counts`: the true unspliced RNA counts (genes x cells)

```{r plot-velocity, fig.width = 4, fig.height = 3.5, out.width = "60%"}
plot_rna_velocity(results, arrow.length = 2)
```

We can inspect the gene-gene correlation using `plot_gene_module_cor_heatmap(results)`:
```{r plot-gene-correlation, fig.width = 8, fig.height = 8}
plot_gene_module_cor_heatmap(results)
```

# Adding Technical Variation and Batch Effect

We can also add the technical variation and batch effect to the true counts.

## Adding technical noise

Simply use the `add_expr_noise` function to add technical noise to the dataset.

```{r technical-noise}
add_expr_noise(
  results,
  # options go here
  alpha_mean = 1e4
)
```

A `counts_obs` field will be added to the `results` object.

This function also accepts a list of options. See the documentation for more details.

- `protocol`: `"umi"` or `"nonUMI"`, whether simulate the UMI protocol.
- `alpha_mean`, `alpha_sd`: Mean and deviation of rate of subsampling of transcripts during capture step.
- `alpha_gene_mean`, `alpha_gene_sd`: `alpha` parameters, but gene-wise.
- `depth_mean`, `depth_sd`: Mean and deviation of sequencing depth.
- `gene_len`: A vector with lengths of all genes.
- `atac.obs.prob`: For each integer count of a particular region for a particular cell, the probability the count will be observed.
- `atac.sd.frac`: The fraction of ATAC-seq data value used as the standard deviation of added normally distrubted noise.
- `randseed`: random seed.

## Adding batch effects

Finally, use the `divide_batches` function to add batch effects.

```{r batch-effects}
divide_batches(
  results,
  nbatch = 2,
  effect = 1
)
```

A `counts_with_batches` field will be added to the `results` object.

The available options are:

- `nbatch`: Number of batches.
- `effect`: The batch effect size.

We can visualize the result with technical noise and batches:

```{r add-expr-noise, fig.width = 4, fig.height = 3.5, out.width = "60%"}
plot_tsne(log2(results$counts_with_batches + 1),
          results$cell_meta$pop,
          legend = 'pop', plot.name = 'RNA Counts Tsne with Batches')
```

# Adjusting Parameters

scMultiSim provides various parameters to control each type of biological effect.
Here, we describe the most important parameters and how they affect the simulation results:

- `num.cifs`, `diff.cif.fraction`
- `cif.mean`, `cif.sigma`
- `discrete.cif`
- `intinsic.noise`

For a complete list of parameters, please check out the [Parameter Guide](https://zhanglabgt.github.io/scMultiSim/articles/options)
page in the documentation.

## The Shiny App

scMultiSim provides a Shiny app to help you generate the options list and visualize the effects of different parameters.
It is highly recommended to use the Shiny app to explore the available parameters.
You can run the app by calling `run_shiny()`.

```{r run-shiny, eval=FALSE}
run_shiny()
```

![Shiny App](https://github.com/ZhangLabGT/scMultiSim/raw/img/img/shiny_app_sc.png)

## Deciding Number of CIFs: `num.cifs`

In scMultiSim, user use `num.cifs` to control the total number of diff-CIF and non-diff-CIFs.
The number of CIFs should be large enough to represent the cell population structure and gene information.
By default, `num.cifs` is set to 50, which is a good starting point for most cases.
However, each gene's base expression is affected by two random diff-CIF entries,
therefore if you have a large number of genes, they may have similar expression patterns, which may not be ideal.
It is recommended to increase `num.cifs` to 50-100 if you have more than 2000 genes.
If you have a small number of genes (less than 1000), you can also decrease `num.cifs` to 20-40.

## Discrete Cell Population: `discrete.cif`

We can also simulate discrete cell population by setting `discrete.cif = TRUE`.
In this case, each tip of the tree will be one cell type,
therefore there will be 5 clusters in the following result.

```{r simulate-discrete, fig.width = 4, fig.height = 3.5, out.width = "60%"}
set.seed(42)

options <- list(
  GRN = GRN_params,
  num.cells = 400,
  num.cifs = 20,
  tree = Phyla5(),
  diff.cif.fraction = 0.8,
  discrete.cif = TRUE
)

results <- sim_true_counts(options)

plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'True RNA Counts Tsne')
```

## Adjusting the Effect of Cell Population: `diff.cif.fraction`

In scMultiSim, the differentiation tree provides explicit control of the cell population.
The effect of the tree can be adjusted by the option `diff.cif.fraction`,
which controls how many CIFs are affected by the cell population.
With a larger `diff.cif.fraction`, the effect of cell population will be larger
and you may see a clearer trajectory or well separated clusters.
With a smaller `diff.cif.fraction`, the resulting RNA counts will be more affected by
other factors, such as the GRN.

Now let's visualize the trajectory with different `diff.cif.fraction` values:

```{r adjust-diff-cif-fraction, fig.width = 4, fig.height = 3.5, out.width = "60%"}
set.seed(42)

options <- list(
  GRN = GRN_params,
  num.cells = 300,
  num.cifs = 20,
  tree = Phyla5(),
  diff.cif.fraction = 0.8
)

results <- sim_true_counts(
        options %>% list_modify(diff.cif.fraction = 0.4))
plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'RNA Counts (diff.cif.fraction = 0.2)')

results <- sim_true_counts(
        options %>% list_modify(diff.cif.fraction = 0.9))
plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'RNA Counts (diff.cif.fraction = 0.8)')
```

## Adjusting the Inherent Cell Heterogeneity: `cif.mean` and `cif.sigma`

The inherent cell heterogeneity is controlled by the non-diff-CIF,
which is sampled from a normal distribution with mean `cif.mean` and standard deviation `cif.sigma`.
Therefore, the larger `cif.sigma` is, the larger the inherent cell heterogeneity is.

Now, let's visualize the effect of `cif.sigma`:

```{r adjust-cif-sigma, fig.width = 4, fig.height = 3.5, out.width = "60%"}
set.seed(42)

options <- list(
  GRN = GRN_params,
  num.cells = 300,
  num.cifs = 20,
  tree = Phyla5(),
  diff.cif.fraction = 0.8,
  cif.sigma = 0.5
)

results <- sim_true_counts(
        options %>% list_modify(cif.sigma = 0.1))
plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'RNA Counts (cif.sigma = 0.1)')

results <- sim_true_counts(
        options %>% list_modify(cif.sigma = 1.0))
plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'RNA Counts (cif.sigma = 1.0)')
```

## Adjusting the Intrinsic Noise: `intinsic.noise`

If we set `do.velocity = FALSE`, scMultiSim will simulate the RNA counts using the Beta-Poisson model,
which is faster but doesn't output RNA velocity.
When using the Beta-Possion model, scMultiSim provides a `intrinsic.noise` parameter to control the
intrinsic noise during the transcription process.
By default, `intrinsic.noise` is set to 1, which means the true counts will be sampled from the Beta-Poisson
model. If we set `intrinsic.noise` to a smaller value like 0.5,
the true counts will be 0.5 * (theoretical mean) + 0.5 * (sampled from the Beta-Poisson model).
_More intrinsic noise will make the encoded effects (e.g. GRN) harder to be inferred._


```{r adjust-intrinsic-noise, fig.width = 4, fig.height = 3.5, out.width = "60%"}
set.seed(42)

options <- list(
  GRN = GRN_params,
  num.cells = 300,
  num.cifs = 20,
  tree = Phyla5(),
  diff.cif.fraction = 0.8,
  intrinsic.noise = 1
)

results <- sim_true_counts(
        options %>% list_modify(intrinsic.noise = 0.5))
plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'RNA Counts (intrinsic.noise = 0.5)')

results <- sim_true_counts(
        options %>% list_modify(intrinsic.noise = 1))
plot_tsne(log2(results$counts + 1),
         results$cell_meta$pop,
         legend = 'pop', plot.name = 'RNA Counts (intrinsic.noise = 1)')
```

## Adjust the effect of chromatin accessibility: `atac.effect`

`atac.effect` Controls the contribution of the chromatin accessibility.
A higher `atac.effect` means the RNA counts are more affected by the ATAC-seq data,
therefore the correlation between the ATAC-seq and RNA-seq data will be higher.

# Simulating Dynamic GRN

First, call the following function to check the usage of dynamic GRN.
```{r help-dynamic-grn}
scmultisim_help("dynamic.GRN")
```

Here we use `Phyla1()` as the differentiation tree to remove the effect of the trajectory. Additionally, we can use `discrete.cif = TRUE` to simulate discrete cell population.
```{r define-options-dynamic-grn}
set.seed(42)

options_ <- list(
  GRN = GRN_params,
  num.cells = 300,
  num.cifs = 20,
  tree = Phyla1(),
  diff.cif.fraction = 0.8,
  do.velocity = FALSE,
  dynamic.GRN = list(
    cell.per.step = 3,
    num.changing.edges = 5,
    weight.mean = 0,
    weight.sd = 4
  )
)

results <- sim_true_counts(options_)
```

`results$cell_specific_grn` is a list containing the gene effects matrix for each cell. Each row is a target and each column is a regulator. The corresponding gene names are displayed as column and row names.
```{r show-cell-specific-grn}
# GRN for cell 1 (first 10 rows)
results$cell_specific_grn[[1]][1:10,]
```

Since we set `cell.per.step = 3`, we expect each adjacent 3 cells share the same GRN:
```{r check-cell-specific-grn}
print(all(results$cell_specific_grn[[1]] == results$cell_specific_grn[[2]]))
print(all(results$cell_specific_grn[[2]] == results$cell_specific_grn[[3]]))
print(all(results$cell_specific_grn[[3]] == results$cell_specific_grn[[4]]))
```

# Session Information

```{r session-info}
sessionInfo()
```