Running a simulation based on SRT data in a single-region • sCCIgen

This tutorial will run a simulation with sCCIgen without the interactive interface based on a SRT data in a single-region.

1 Load R package

library(sCCIgen)

2 Load and clean sample data

Download sample data from https://github.com/songxiaoyu/sCCIgen_data/tree/main/input_data

load("SeqFishPlusCortex_2025_expr.Rdata")
load("SeqFishPlusCortex_2025_spatial.Rdata")

dim(expr)
# [1] 2500  511
dim(spatial)
# [1] 511   3
expr[1:3, 1:3]
#       ExcitatoryNeuron ExcitatoryNeuron ExcitatoryNeuron
# Aatf                 0                0                0
# Abat                 1                0                0
# Abhd2                2                2                0
spatial[1:3,]
#               anno X_final  Y_final
# 1 ExcitatoryNeuron 1632.02 -1305.70
# 2 ExcitatoryNeuron 1589.47  -669.51
# 3 ExcitatoryNeuron 1539.89 -1185.90
anno <- colnames(expr)

3 Analysis of the existing data to provide insights into the parameters of the simulation

Users can split the sCCIgen simulation into (1) model fitting and (2) simulation using fitted model and user-provided parameters steps to expedite the simulations.

It is especially helpful if the number of genes and/or cells are very large and users want to run simulation for more than once. By splitting the simulation into these two steps, users can estimate model parameters only once and save the results for multiple use.

3.1 Task 1: Estimate model parameters from the SRT data for simulation

This is part is to fit the expression data. When sim_method=“copula”, it will fit both the gene marginal distribution and gene-gene correlation. When sim_method=“ind”, it will only fit the gene marginal distribution.

Note: If the number of genes or cells are large, model fitting may take some time. It is suggested to select a reasonable sample size (e.g. <1500 cells per cell type) before the model fitting, as more cells may not be needed improve the estimation. The user interface automatically uses a maximum of 2500 cells per cell type in model fitting. In the meantime, if some genes are extremely zero-inflated, narrowing the simulation to reasonably variable genes is an option.

# model fitting w/o gene-gene correlation (fast but less accurate estimation)
ModelEst <- Est_ModelPara(expr = expr, 
                          anno = anno, 
                          sim_method = "ind",
                          ncores = 10)

saveRDS(ModelEst, file = "SeqFishPlusCortex_fit_wo_cor.RDS")

# model fitting with gene-gene correlation
ModelEst <- Est_ModelPara(expr = expr, 
                          anno = anno, 
                          sim_method = "copula", 
                          ncores = 10)

saveRDS(ModelEst, file = "SeqFishPlusCortex_fit_w_cor.RDS")

3.2 Task 2: Estimate CCIs in the input data following the existing Giotto pipeline

library(Giotto)
# Preprocess data with Giotto pipeline
## Three types of networks are built in to determine cell neighborhood: Delaunay, kNN, and distance-based. Here we use distance-based network as an example. 

db <- preprocessGiotto(expr_data = expr, 
                       spatial_data = spatial, 
                       run_hvg = TRUE,
                       run_Dist_network = TRUE,
                       dis.cut=150)
db
# Estimate cell-cell attraction and inhibition patterns, and save in pre-defined folder
cellProximityTable(gobject = db, 
                   abs_enrichm = 0.3, 
                   p_adj = 0.05,
                   spatial_network_name = "distance_based_network",
                   save.unfiltered=TRUE,
                   output_file = "est_CCI_dist_dist.csv",
                   seed=123)

# Estimate gene expressions of cells impacted by their neighbors, and save in pre-defined folder
ExprDistanceTable(gobject = db, 
                  in_hvg = TRUE, 
                  region_specific = FALSE, 
                  spatial_network_name = "distance_based_network",
                  abs_log2fc_ICG = 0.3, 
                  p_adj = 0.05,
                  output_file = "est_CCI_dist_expr.csv",
                  seed=123)

# Estimate gene expressions of cells impacted by gene expressions of neighboring cells, narrow to
#  known pairs, such as ligand and receptor pairs, and save in pre-defined folder
ExprExprTable(gobject = db, 
              database = "mouse", 
              region_specific = FALSE, 
              spatial_network_name = "distance_based_network",
              p_adj = 0.05, 
              abs_log2fc_LR = 0.2,
              output_file = "est_CCI_expr_expr.csv",
              seed=123)

4 Create a parameter file

Users need to create a parameter file. The sample parameter file for SRT based simulation is here for downloading and filling in to perform simulations.

5 Perform the entire simulation and save the results.

Assuming you already have a parameter file, you can run the entire simulation using codes like this:

# Simulate default data - using existing cells but simulate expression with ground truth
input <- "SeqFishPlus_default_param.yml"
model_param_path="SeqFishPlusCortex_fit_wo_cor.RDS"
ParaSimulation(input = input, ModelFitFile = model_param_path)



# Simulate new cells of the same number (n=511) and expression of genes with gene-gene correlation expression, including three types of estimated CCIs.
input <- "SeqFishPlus_n511_cor_CCIs_param.yml"
ParaSimulation(input = input)