dicepro - Simulated Data Workflow

dicepro Team

2026-06-24

Note: All code chunks have eval = FALSE and are shown for illustration only. To run them interactively:

library(dicepro)
# copy-paste the chunks below into your R session

1 Overview

This vignette demonstrates complete dicepro deconvolution workflows on fully synthetic data, enabling end-to-end reproducibility without requiring access to proprietary data-sets.

Two simulation strategies are available and covered here:

Function Reference matrix Cell-type names Best used for
simulation() Synthetic (generated internally) Generic (CellType_k) Flexible benchmarking with custom parameters
simulation_bluecode() BlueCode (bundled real reference) Real biological names Realistic benchmarking on human tissue cell types

Both functions return the same standardized triplet – $W, $p, $B – and are fully interchangeable as input to dicepro().

The pipeline proceeds in four stages:

  1. Simulate a reference signature matrix and ground-truth cell-type proportions.
  2. Generate a noisy bulk RNA-seq matrix from the simulated mixture.
  3. Deconvolve using the main dicepro() function with automated hyper-parameter optimization.
  4. Evaluate the recovered proportions against the ground truth.

2 Strategy 1 – Fully Synthetic Simulation

simulation() generates a self-consistent triplet from scratch: a synthetic reference matrix $W, ground-truth proportion matrix $p, and noisy bulk matrix $B, guaranteeing column-name alignment across all three outputs.

2.1 Data Generation

library(dicepro)
set.seed(2101L)

sim_data <- simulation(
  loi        = "gauss",
  scenario   = "hierarchical",
  nSample    = 30,
  nGenes     = 200,
  nCellsType = 10,
  sigma_bio  = 0.07,
  sigma_tech = 0.07,
  seed       = 2101L
)

cat("Reference  :", dim(sim_data$W), "\n")
cat("Proportions:", dim(sim_data$p), "\n")
cat("Bulk       :", dim(sim_data$B), "\n")
cat("Row sums   :", round(range(rowSums(sim_data$p)), 4), "\n")

2.2 Noise Model Sanity Check

bulk_clean <- as.matrix(sim_data$W) %*% t(as.matrix(sim_data$p))

plot(
  bulk_clean[seq_len(500)],
  as.matrix(sim_data$B)[seq_len(500)],
  xlab = "Clean bulk (first 500 entries)",
  ylab = "Noisy bulk",
  pch  = 19, cex = 0.4, col = "#2c7bb6",
  main = "Noise model: clean vs noisy bulk"
)
abline(0, 1, col = "firebrick", lwd = 1.5)

2.3 Deconvolution

out <- dicepro(
  reference             = as.matrix(sim_data$W)[, -c(1,5,10)],
  bulk                  = as.matrix(sim_data$B),
  methodDeconv          = "FARDEEP",
  W_prime               = 0,
  bulkName              = "SimBulk",
  refName               = "SimRef",
  hp_max_evals          = 500,
  algo_select           = "random",
  output_path           = tempdir(),
  hspaceTechniqueChoose = "gamma_dominant",
  normalize             = FALSE
)

2.4 Results

2.4.1 Optimal Hyperparameters

cat("Best lambda :", out$hyperparameters$lambda, "\n")
cat("Best gamma  :", out$hyperparameters$gamma,  "\n")
cat("Loss        :", out$metrics$loss,           "\n")
cat("Constraint  :", out$metrics$constraint,     "\n")

2.4.2 Hyperparameter Optimisation Report

out$plot_hyperopt

2.4.3 Pareto Frontier

out$plot

2.4.4 Recovered vs True Proportions

true_prop <- as.matrix(sim_data$p)
pred_prop <- as.matrix(out$H)

common_ct   <- intersect(colnames(pred_prop), colnames(true_prop))
true_common <- true_prop[, common_ct, drop = FALSE]
pred_common <- pred_prop[, common_ct, drop = FALSE]

r_overall <- cor(as.vector(true_common), as.vector(pred_common))
cat(sprintf("Overall Pearson r: %.3f\n", r_overall))

plot(
  as.vector(true_common),
  as.vector(pred_common),
  xlab = "True proportions",
  ylab = "Predicted proportions",
  pch  = 19, cex = 0.5, col = "#2c7bb6aa",
  main = sprintf("True vs Predicted  (r = %.3f)", r_overall)
)
abline(0, 1, col = "firebrick", lwd = 1.5)

2.4.5 Per-Cell-Type Correlation

ct_cors <- vapply(common_ct, function(ct) {
  cor(true_common[, ct], pred_common[, ct])
}, numeric(1L))

par(mar = c(5, 10, 3, 1))
barplot(
  sort(ct_cors),
  horiz = TRUE, las = 1,
  col   = ifelse(sort(ct_cors) > 0.7, "#2c7bb6", "#d7191c"),
  xlab  = "Pearson r",
  main  = "Per-cell-type correlation (synthetic)",
  xlim  = c(-0.2, 1)
)
abline(v = 0.7, lty = 2, col = "gray40")

2.4.6 Quantitative Performance Metrics

perf <- makeTable1Tool(pred_mat = pred_common, obs_mat = true_common)
knitr::kable(perf$Perf, digits = 3,
             caption = "Performance metrics -- fully synthetic data")

3 Strategy 2 – BlueCode-Based Simulation

simulation_bluecode() uses the BlueCode reference matrix bundled with the package – a curated signature matrix covering 34 human cell types across five tissue compartments – as the fixed reference $W. Proportions are simulated via a two-level Dirichlet hierarchy that reflects the biological organisation of human tissues.

This strategy is recommended when benchmarking should be grounded in a biologically realistic reference, while keeping full control over the ground-truth proportions.

3.1 Compartment Structure

The hierarchy encodes five major tissue compartments with their constituent cell types:

Compartment Cell types (n) Default α
Immune B naive/memory, T CD4/CD8, NK, Monocyte, Macrophage, Neutrophil 4.0
Stromal Fibroblasts (×5), MSC, Chondrocyte, Osteoblast 2.5
Endothelial Large vessel, microvascular (×2) 1.8
Epithelial Mammary, renal, respiratory, keratinocyte, melanocyte 1.8
Muscle Smooth muscle (×7), cardiac myocyte, myometrial 1.5

3.2 Data Generation

library(dicepro)

sim_bc <- simulation_bluecode(
  nSample    = 30,
  sigma_bio  = 0.15,
  sigma_tech = 0.02,
  seed       = 2101L
)

cat("Reference  :", dim(sim_bc$W), "\n")  # nGenes x 34
cat("Proportions:", dim(sim_bc$p), "\n")  # 30 x 34
cat("Bulk       :", dim(sim_bc$B), "\n")  # nGenes x 30
cat("Row sums   :", round(range(rowSums(sim_bc$p)), 4), "\n")

# Real cell-type names from BlueCode
head(colnames(sim_bc$p))

3.3 Proportion Distribution by Compartment

The hierarchical model induces realistic between-compartment variation. Visualising the compartment-level totals confirms the expected dominance of the Immune compartment.

# Cell-type to compartment mapping (mirrors .bluecode_cell_groups)
compartment_map <- c(
  rep("Immune",      9),
  rep("Stromal",     8),
  rep("Endothelial", 3),
  rep("Epithelial",  5),
  rep("Muscle",      9)
)
names(compartment_map) <- colnames(sim_bc$p)

# Aggregate proportions by compartment for each sample
comp_props <- t(apply(sim_bc$p, 1, function(row) {
  tapply(row, compartment_map[names(row)], sum)
}))

boxplot(
  comp_props,
  col  = c("#4393c3", "#74c476", "#fd8d3c", "#9ecae1", "#fb6a4a"),
  ylab = "Compartment proportion",
  main = "Simulated compartment proportions (BlueCode)",
  las  = 2
)

3.4 Noise Model Sanity Check

bulk_clean_bc <- as.matrix(sim_bc$W) %*% t(as.matrix(sim_bc$p))

plot(
  bulk_clean_bc[seq_len(500)],
  as.matrix(sim_bc$B)[seq_len(500)],
  xlab = "Clean bulk (first 500 entries)",
  ylab = "Noisy bulk",
  pch  = 19, cex = 0.4, col = "#74c476",
  main = "Noise model: clean vs noisy bulk (BlueCode)"
)
abline(0, 1, col = "firebrick", lwd = 1.5)

3.5 Deconvolution

out_bc <- dicepro(
  reference             = as.matrix(sim_bc$W),
  bulk                  = as.matrix(sim_bc$B),
  methodDeconv          = "FARDEEP",
  W_prime               = 0,
  bulkName              = "BlueCodeBulk",
  refName               = "BlueCode",
  hp_max_evals          = 100,
  algo_select           = "random",
  output_path           = tempdir(),
  hspaceTechniqueChoose = "gamma_dominant",
  normalize             = FALSE
)

3.6 Results

3.6.1 Optimal Hyperparameters

cat("Best lambda :", out_bc$hyperparameters$lambda, "\n")
cat("Best gamma  :", out_bc$hyperparameters$gamma,  "\n")
cat("Loss        :", out_bc$metrics$loss,            "\n")
cat("Constraint  :", out_bc$metrics$constraint,      "\n")

3.6.2 Hyperparameter Optimisation Report

out_bc$plot_hyperopt

3.6.3 Pareto Frontier

out_bc$plot

3.6.4 Recovered vs True Proportions

true_prop_bc <- as.matrix(sim_bc$p)
pred_prop_bc <- as.matrix(out_bc$H)

common_ct_bc   <- intersect(colnames(pred_prop_bc), colnames(true_prop_bc))
true_common_bc <- true_prop_bc[, common_ct_bc, drop = FALSE]
pred_common_bc <- pred_prop_bc[, common_ct_bc, drop = FALSE]

r_overall_bc <- cor(as.vector(true_common_bc), as.vector(pred_common_bc))
cat(sprintf("Overall Pearson r: %.3f\n", r_overall_bc))

plot(
  as.vector(true_common_bc),
  as.vector(pred_common_bc),
  xlab = "True proportions",
  ylab = "Predicted proportions",
  pch  = 19, cex = 0.5, col = "#74c47699",
  main = sprintf("True vs Predicted -- BlueCode  (r = %.3f)", r_overall_bc)
)
abline(0, 1, col = "firebrick", lwd = 1.5)

3.6.5 Per-Cell-Type Correlation

ct_cors_bc <- vapply(common_ct_bc, function(ct) {
  cor(true_common_bc[, ct], pred_common_bc[, ct])
}, numeric(1L))

par(mar = c(5, 14, 3, 1))
barplot(
  sort(ct_cors_bc),
  horiz = TRUE, las = 1,
  col   = ifelse(sort(ct_cors_bc) > 0.7, "#2c7bb6", "#d7191c"),
  xlab  = "Pearson r",
  main  = "Per-cell-type correlation (BlueCode)",
  xlim  = c(-0.2, 1)
)
abline(v = 0.7, lty = 2, col = "gray40")

3.6.6 Quantitative Performance Metrics

perf_bc <- makeTable1Tool(pred_mat = pred_common_bc, obs_mat = true_common_bc)
knitr::kable(perf_bc$Perf, digits = 3,
             caption = "Performance metrics -- BlueCode simulation")

4 Comparing Both Strategies

When both strategies are run under the same conditions, their overall correlation scores can be compared directly.

cat(sprintf(
  "Overall Pearson r\n  Synthetic  : %.3f\n  BlueCode   : %.3f\n",
  r_overall, r_overall_bc
))

These two strategies are complementary: simulation() offers maximum flexibility for testing deconvolution under stress conditions, with varying numbers of genes, cell types, or noise levels; simulation_bluecode() anchors the benchmark to a biologically grounded reference, making the results easier to interpret in the context of real-world data from human tissues.

5 Session Info

sessionInfo()