Simulating Single-Cell Sequencing Data at the Clonal Resolution

Introduction

This vignette demonstrates how to use the dynamics module to simulate single-cell dynamics, visualize population history, and generate lineage trees. The module is highly customizable and can be adapted to various biological scenarios. Key features include the ability to simulate cell birth, death, and transition events, track mutations, and analyze population dynamics over time.

Dynamics Parameters

The simulation relies on several key parameters to model single-cell dynamics:

Initial Population: A named vector specifying the starting number of cells in each clone.
Birth and Death Rates: Intrinsic rates governing cell division and death for each clone.
Transition Rates: Rates at which cells transition from one clone to another, defined in a data frame with parent, child, and rate columns.
Clone Capacity: The maximum population size each clone can reach, simulating resource limitations.
Chromosome Lengths: A list specifying the lengths of chromosomes for each clone, used to model mutations.
Mutation Rate: The probability of a mutation occurring during cell division.

These parameters allow users to fine-tune the simulation to reflect specific biological systems or experimental conditions.

Simulation Scenario

In this vignette, we demonstrate the use of Scenario 2 from our paper (Jin et al., 2024, bioRxiv), which models the emergence of a single dominant clone within a population. This scenario is characterized by:

A high birth rate for the dominant clone (e.g., 1.0) compared to other subclones (e.g., 0.5).
A low death rate across all clones, ensuring sustained population growth.
A limited carrying capacity for each clone, simulating resource constraints in the environment.

To evaluate the robustness of the simulation, we ran three replicates of this scenario. Despite the inherent stochasticity of the underlying Markov process, all replicates exhibited consistent population dynamics, with the dominant clone achieving prominence while other subclones remained in the minority. This reproducibility highlights the reliability of our simulation framework for modeling evolutionary dynamics.

Load Data and Set Parameters

Load simulation data and define parameters for the simulation.

# Load libraries
library(cancerSimCraft)
library(igraph)
library(gridExtra)

# Load simulation data
load(file = "./tutorials_data/advanced_4_simulation_setting.RData")
load(file = "./large_tutorials_data/advanced_6_clonal_level_cancer_genome_simulation.RData")

# Set initial parameters
initial_birth_rates <- c(0.05, 0.5, 0.5, 0.5, 0.5, 0.5, 1)
initial_death_rates <- 0.3 * c(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01)
initial_population <- c(1, 0, 0, 0, 0, 0, 0)
clone_capacity <- c(500, 2000, 5000, 5000, 5000, 5000, 5000)

# Name parameters
names(initial_population) <- names(initial_birth_rates) <- names(initial_death_rates) <- names(clone_capacity) <- V(chisel_tree)$name

# Define transition rates
initial_edge_transition_rates <- data.frame(as_edgelist(chisel_tree), rate = rep(0.0001, 6))
colnames(initial_edge_transition_rates) <- c("parent", "child", "rate")

Run Simulations

Run multiple simulations and store the results.

# Define the number of simulation sets
num_simulation_sets <- 3
max_steps <- 20000

# Create a list to store simulation results
param2_dynamics_sim <- list()

# Run simulations
for (sim_set in 1:num_simulation_sets) {
  dynamics_sim <- simulate_sc_dynamics(
    initial_population = initial_population,
    max_steps = max_steps,
    intrinsic_birth_rates = initial_birth_rates,
    intrinsic_death_rates = initial_death_rates,
    edge_transition_rates = initial_edge_transition_rates,
    clone_capacity = clone_capacity,
    chr_lengths = sim_updates$all_node_chr_lengths,
    mutation_rate = 1
  )
  param2_dynamics_sim[[sim_set]] <- dynamics_sim
}

Visualize Population History

Generate and save population history plots for each simulation.

# Define clone colors
goodColor2 <-  c("black","#DF536B","#61D04F","#2297E6","#28E2E5","#CD0BBC","#F5C710","gray62")
clone_colors <- goodColor2[1:length(clone_names)]
names(clone_colors) <- clone_names

# Generate population history plots
plot_list <- list()
for (sim_set in 1:num_simulation_sets) {
  title <- paste0("Scenario2, Rep", sim_set)
  plot_list[[sim_set]] <- plot_population_history(param2_dynamics_sim[[sim_set]], clone_colors, title = title)
  print(plot_list[[sim_set]])
}

Visualize Lineage Trees

Generate and save lineage tree plots for each simulation.

# Generate lineage tree plots
tree_plot_list <- list()
for (sim_set in 1:num_simulation_sets) {
  title <- paste0("Scenario2, Rep", sim_set)
  tree_plot_list[[sim_set]] <- create_sc_tree_plot(
    dynamic_sim_ob = param2_dynamics_sim[[sim_set]],
    clone_colors = clone_colors,
    title = title
  )
  print(tree_plot_list[[sim_set]])
}

Save Dynamics Data

Save the final dynamics data for future analysis.

save.image(file = "./tutorials_data/advanced_8_dynamics_simulation.RData")