SimiCviz: Visualization and Analysis of Gene Regulatory Networks
3 August 2026
Contents
1 Introduction
SimiCviz provides various utilities for gene regulatory network (GRN) analysis for single-cell RNA-seq data. Originally developed to visualize the phenotype-specific GRN output from SimiCPipeline for R users, it also includes tools for:
- Load and organize GRN outputs from multiple sources (SimiCPipeline, SCENIC, Pando, etc.)
- Calculate cell activity scores (AS) using expression data and GRN weights.
- Visualize regulatory networks, distribution of TF/regulon activities, and dissimilarity across conditions.
The package is interoperable and agnostic to the GRN inference methods, accepting different GRN outputs.
1.1 Relationship to Similar Packages
Packages such as SCENIC and Pando focus on reconstructing or modelling
gene regulatory networks from single-cell data, often producing regulons,
TF-target links, activity estimates, or statistical evidence for regulatory
interactions. SimiCviz is complementary to these tools. It does not replace
their inference workflows or duplicate their visualization implementations;
instead, it provides a common downstream interface for loading GRN outputs,
organizing TF-target weights and cell annotations, computing compatible activity
scores when expression data are available, and generating comparative
visualizations across phenotypes or cell groups.
For SimiCPipeline outputs, SimiCviz provides direct readers for the
standard directory structure and pickle/CSV files. For other methods, including
SCENIC, Pando, or custom GRN workflows, users can import a long-format table
with transcription factors, targets, weights, and optional quality metrics such
as adjusted p-values or model-fit statistics.
2 Installation
The development version can be installed from GitHub:
if (!requireNamespace("remotes", quietly = TRUE)) {
install.packages("remotes")
}
remotes::install_github("ML4BM-Lab/SimiCviz")Some import helpers use reticulate to read Python-backed formats such as
pickle or H5AD files. When those functions are called, SimiCviz requests a
Python environment through reticulate::py_require() with Python >= 3.8 and
the numpy, pandas, and anndata Python packages. These Python dependencies
are only needed for those import paths; CSV and R-native workflows do not
require them.
The package will be available from Bioconductor upon approval:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("SimiCviz")3 Loading GRN Data
3.1 Example Dataset
The examples in this vignette use a reduced dataset derived from the
SimiC-Suite case study.
The original dataset is a single-cell RNA-seq dataset of multiple myeloma (MM) patients, including healthy donors and patients at different disease stages (Boiarsky et al. 2022).
The raw data is available from the GEO repository (GSE193531).
We followed the workflow described in the tutorial and then generated a smaller subset for Bioconductor vignette examples, keeping representative
SimiCPipeline outputs for the NBM, SMM, and MM disease-stage labels while
limiting the package size and vignette runtime. Exact code used is included in inst/scripts/generate_example.ipynb and generate_simicfull.R.
3.2 Data Format Requirements
SimiCviz easily handles SimiCPipeline outputs but also accepts GRN weights in a flexible long format. The input data.frame must include the following columns:
tf: Transcription factor nametarget: Target gene name
weight: Edge weight/coefficientlabel(optional): Phenotype/condition identifier for label-specific networks
Optional quality metrics for filtering weights:
pvalue/adj_p_val: Statistical significance (for methods like SCENIC, Pando)r_squared/adj_r_squared: Goodness of fit (for methods like SimiC)
3.3 SimiCPipeline Outputs
SimiCPipeline is a regularized regression method for inferring condition-specific GRNs. It generates:
.picklefiles with TF-target weights per phenotype/label and adjusted R² values for goodness of fit of the target expression model..csvfile with a TF activity score matrix (cells × TFs).picklefiles with per-label matrices (cells × TFs) [Legacy SimiC v1 - Deprecated]
These outputs are organized in a directory structure like this:
Project/
├── inputFiles/
│ ├── TF_list.csv
│ ├── expression_matrix.pickle
│ └── phenotype_annotation.txt
└── outputSimic/
├── figures/
└── matrices/
└── example1/
├── example1_L1_0.1_L2_0.01_simic_matrices.pickle
├── example1_L1_0.1_L2_0.01_simic_matrices_filtered_BIC.pickle
├── example1_L1_0.1_L2_0.01_wAUC_matrices_filtered_BIC.pickle
└── example1_L1_0.1_L2_0.01_wAUC_matrices_filtered_BIC_collected.csvHere we show how to easily load a SimiCPipeline run. You only need the directory, experiment name, and hyperparameters since the directory output is standardized. SimiCviz will automatically find the files and generate the SimiCvizExperiment.
library(SimiCviz)
# Load entire SimiCPipeline run automatically
simic_full <- load_SimiCPipeline(
project_dir = "path/to/simic_run",
run_name = "example1",
lambda1 = "0.01",
lambda2 = "0.001"
)
# Set display names and colors for visualization (Part 3)
simic_full <- setLabelNames(
simic_full,
label_names = c('NBM', 'SMM', 'MM'),
colors = c("#3B7EA1", "#E66101", "#B2182B")
)simic_full
## An object of class SimiCvizExperiment
## 3 label(s), 10 TF(s), 150 target(s)
## Weights: 3 matrices [0: 10 x 150, 1: 10 x 150, 2: 10 x 150]
## AUC: collected (2250 cells x 10 TFs)
## Cell labels: 2250 cells across 3 label(s) [0, 1, 2]
## Label names: 0 = NBM, 1 = SMM, 2 = MM
## Colors: 0 = #3B7EA1, 1 = #E66101, 2 = #B2182B
## TFs: MEF2D, E2F4, SATB1, IRF1, ATF6, JUN, ...
## Targets: RPLP1, RPL36, RPL39, IGHM, EEF1A1, RPS8, ...
## Meta keys: adjusted_r_squaredIf you run the full SimiCPipeline tutorial you can skip Part 2 and go straight to Part 3 for visualizations.
3.4 Manual Loading
If you need more flexibility in the workflow, or have precomputed outputs using SimiC v1, you can construct SimiCvizExperiment or AUCProcessor objects by loading the required files separately.
3.4.1 Weights from pickle
# Load weights from pickle
weights_file <- system.file("extdata",
file.path("outputSimic/example_simic_weights.pickle"),
package = "SimiCviz")
simic_weights <- read_weights_pickle(weights_file)
simic_weights[[1]][, 1:6]
## RPLP1 RPL36 RPL39 IGHM
## MEF2D -0.7501612 0.0000000 0.0000000 0.0000000
## E2F4 0.0000000 0.0000000 0.0000000 0.9553137
## SATB1 0.0000000 0.0000000 0.0000000 0.0000000
## IRF1 -0.8536968 0.0000000 -0.8289669 1.1979259
## ATF6 -2.6133661 -2.2118857 -2.0251983 0.0000000
## JUN 0.0000000 0.0000000 0.0000000 -0.8878323
## JUND 0.0000000 0.0000000 0.0000000 0.0000000
## POU2F2 0.0000000 0.7736768 0.0000000 0.0000000
## KLF3 0.0000000 0.0000000 0.0000000 0.0000000
## BCL11A 0.7300538 0.0000000 0.0000000 0.9884850
## EEF1A1 RPS8
## MEF2D 0.0000000 0.000000
## E2F4 0.0000000 0.000000
## SATB1 0.0000000 0.000000
## IRF1 0.0000000 -1.112810
## ATF6 -1.8176760 -2.430652
## JUN 0.0000000 0.000000
## JUND 0.0000000 0.000000
## POU2F2 0.9167308 0.000000
## KLF3 0.0000000 0.000000
## BCL11A 0.0000000 0.0000003.4.2 Weights from other methods
Most single-cell GRN inference methods produce table outputs with TF-target weights/estimates and p-values.
3.4.2.1 Example: Generic CSV Format
# Load GRN weights (method agnostic)
weight_path <- system.file("extdata", "example_weights.csv",
package = "SimiCviz")
# Read as data.frame
weights_df <- read_weights_csv(weight_path)
head(weights_df)
## tf target weight label
## 1 MEF2D RPLP1 -0.7501612 0
## 2 E2F4 RPLP1 0.0000000 0
## 3 SATB1 RPLP1 0.0000000 0
## 4 IRF1 RPLP1 -0.8536968 0
## 5 ATF6 RPLP1 -2.6133661 0
## 6 JUN RPLP1 0.0000000 0
# If your method uses different column names you need to rename themExpected columns in CSV:
# Minimal
tf, target, weight
# If only one GRN a column name "label" will be all values in label should be 0
# With quality metrics
tf, target, weight, pvalue, adj_p_val
# With phenotype labels
tf, target, weight, label, pvalue, adj_p_val3.4.3 Cell Labels / Phenotype Annotations
The cell_labels map each cell ID in your expression matrix or activity scores matrix to a condition/phenotype. This is required for activity score calculation and AUC visualizations.
# Load from CSV (recommended format: columns 'cell', 'label')
cell_labels_path <- system.file("extdata",
file.path("inputFiles", "disease_stage_annotation.csv"),
package = "SimiCviz")
cell_labels <- load_cell_labels(cell_labels_path, header = TRUE, sep = ",")
## Cell labels file contains extra columns
## Cell / category / label
head(cell_labels)
## cell category label
## 1 AGGGTGATCTGAGGGA-1-NBM-10.138P NBM 0
## 2 CCTAGCTTCTCCAACC-1-NBM-1.138P NBM 0
## 3 ATTACTCTCGTGGTCG-1-NBM-10.138P NBM 0
## 4 TAAGAGAAGCCGCCTA-1-NBM-8.138P NBM 0
## 5 GATCGATTCAGAGGTG-1-NBM-8.138P NBM 0
## 6 ACATCAGGTCGCGGTT-1-NBM-10.138P NBM 0The function load_cell_labels can take alternative formats/files:
# From vector (will generate cell_1, cell_2, ... names)
cell_labels <- c(0, 0, 1, 1, 2, 2) # Must match order of AUC rows
# From named vector
cell_labels <- c(cell_A = 0, cell_B = 0, cell_C = 1)
# From data.frame
cell_labels <- data.frame(cell = c("cell_A", "cell_B"), label = c(0, 1))3.4.4 Expression Matrix
Only required for computing activity scores. Can be:
- CSV file
- Pickle file: numpy array or pandas DataFrame
- H5AD file: AnnData object
- RDS file: Seurat or SingleCellExperiment object
- R object: matrix or data.frame
Important: Final expression matrix must be in genes × cells format (genes as rows, cells as columns).
# Load from multiple formats
expression_mat_path <- system.file("extdata",
file.path("inputFiles", "example_expression.pickle"),
package = "SimiCviz")
# Will auto-detect format and load
expression_mat <- load_expression_matrix(expression_mat_path)
## Converting expression matrix to sparse format...
print(class(expression_mat))
## [1] "dgCMatrix"
## attr(,"package")
## [1] "Matrix"
print(dim(expression_mat))
## [1] 2250 1603.4.5 Create a SimiCvizExperiment Object
The SimiCvizExperiment container organizes weights, activity scores, and metadata.
# From .pickle SimiC files
# Extract Adjusted R squared from SimiC outputs
out <- read_pickle(weights_file)
adjusted_r_squared <- out$adjusted_r_squared
viz_obj_simic <- SimiCvizExperiment(
weights = simic_weights,
auc = NULL, # Will compute this in the next section later but can be loaded as well
cell_labels = cell_labels,
label_names = c('NBM', 'SMM', 'MM'),
colors = c("#3B7EA1", "#E66101", "#B2182B"),
meta=list(adjusted_r_squared=adjusted_r_squared))
## `weights` provided as a list.
## Cell labels file contains extra columns
## cell / category / label
viz_obj_simic
## An object of class SimiCvizExperiment
## 3 label(s), 10 TF(s), 150 target(s)
## Weights: 3 matrices [0: 10 x 150, 1: 10 x 150, 2: 10 x 150]
## AUC collected: none
## Cell labels: 2250 cells across 3 label(s) [0, 1, 2]
## Label names: 0 = NBM, 1 = SMM, 2 = MM
## Colors: 0 = #3B7EA1, 1 = #E66101, 2 = #B2182B
## TFs: MEF2D, E2F4, SATB1, IRF1, ATF6, JUN, ...
## Targets: RPLP1, RPL36, RPL39, IGHM, EEF1A1, RPS8, ...
## Meta keys: adjusted_r_squared4 Computing Activity Scores (AUC)
In many biological contexts, we aim to quantify changes in GRN activity across cell populations. This is achieved by computing an activity score for each TF in each cell, as described in the original SimiC article by Peng et al published in Communications Biology (2022). Here, we provide wrapper functions to apply this approach to outputs from SimiCPipeline or other GRN inference methods.
4.1 Quick Start
Compute activity scores from TF-Target weights and expression matrix using SimiCvizExperiment
# From SimiC input files
# adj_r2_threshold: For SimiCPipeline style outputs it will look in metadata of viz_obj_simic, for "simic@meta$adjusted_r_squared" and filter out targets with lower R²)
# n_cores: Number of workers if backend = 'multissession' or
# Number of cores if backend = 'multicore'
viz_obj_simic <- calculate_activity_scores(
viz_obj_simic,
expression = expression_mat_path,
adj_r2_threshold = 0.7, # For SimiC style outputs
sort_by="expression", # Rank targets by expression or weight
select_top_k = NULL, # Use all targets (or limit to top K)
percent_of_target = 1.0,
n_cores = 2,
backend = "multicore",
verbose = TRUE
)
## Initializing AUCProcessor...
## `weights` provided as a list, converting to list data.frame
## Filtering weights by adjusted R2 threshold...
## Converting expression matrix to sparse format...
## Cell labels file contains extra columns
## cell / category / label
## Transposing expression matrix to
## genes x cells format based on rownames matching cell IDs
## Computing activity scores...
## Starting AUC computation...
## Sorting by: expression
## Targets: 100% of available
## Labels: 3 | Cells: 2250
## Backend requested: multicore | backend used: multicore |
## workers: 2
##
|
| | 0%
|
|====== | 12%
|
|============ | 25%
|
|=================== | 38%
|
|========================= | 50%
|
|=============================== | 62%
|
|====================================== | 75%
|
|============================================ | 88%
|
|==================================================| 100%
## AUC computation completed in 5.37 seconds!
## Result: 2250 cells x 10 TFs
## Activity scores computed and added to simic object
# Access computed scores
auc_scores <- viz_obj_simic@auc$collected
head(auc_scores[, 1:5])
## MEF2D E2F4 SATB1
## AGGGTGATCTGAGGGA-1-NBM-10.138P 0.3390212 0.2978777 NaN
## CCTAGCTTCTCCAACC-1-NBM-1.138P 0.3672169 0.2732215 NaN
## ATTACTCTCGTGGTCG-1-NBM-10.138P 0.3204024 0.2985785 NaN
## TAAGAGAAGCCGCCTA-1-NBM-8.138P 0.3382259 0.2739050 NaN
## GATCGATTCAGAGGTG-1-NBM-8.138P 0.3529123 0.2737598 NaN
## ACATCAGGTCGCGGTT-1-NBM-10.138P 0.3330106 0.2955335 NaN
## IRF1 ATF6
## AGGGTGATCTGAGGGA-1-NBM-10.138P 0.3703713 0.3460459
## CCTAGCTTCTCCAACC-1-NBM-1.138P 0.3649993 0.3376641
## ATTACTCTCGTGGTCG-1-NBM-10.138P 0.3696492 0.3350991
## TAAGAGAAGCCGCCTA-1-NBM-8.138P 0.3536551 0.3438787
## GATCGATTCAGAGGTG-1-NBM-8.138P 0.3477387 0.3518685
## ACATCAGGTCGCGGTT-1-NBM-10.138P 0.3646941 0.33934184.2 Advanced: AUCProcessor Workflow
For more control over the computation pipeline, you can use AUCProcessor class directly.
If your method outputs a quality metric for the inferred GRN, you can filter out unreliable TF-target relationships with the arguments:
qc_type: “adj_r2” or name of the column inweightsto use (pval, padj, etc.). If ‘qc_type’ == “adj_r2” then cells aboveqc_thresholdwill be kept (e.g. R² > 0.7). If ‘qc_type’ != “adj_r2” (we assume is a pvalue) then cells belowqc_thresholdwill be kept (e.g. adj_p_val < 0.05).qc_threshold: desired threshold as described earlier.
Tip: If your method output is different you should modify your input column names accordingly.
# Initialize processor with weights and expression
AS_processor <- AUCProcessor(
weights = weights_df,
expression = expression_mat_path, # or matrix/data.frame
cell_labels = cell_labels,
qc_type = "adj_p_val",
qc_threshold = 0.05, # Filter targets above this threshold
n_cores = 2,
backend = "multisession"
)
## `weights` provided as a data.frame
## No adj_p_val column found in weights data.frame;
## skipping filtering.
## Converting expression matrix to sparse format...
## Cell labels file contains extra columns
## cell / category / label
## Transposing expression matrix to
## genes x cells format based on rownames matching cell IDs
# Compute with custom parameters
AS_processor <- compute_auc(
AS_processor,
sort_by = "expression",
select_top_k = NULL, # Use all targets (or limit to top K)
percent_of_target = 1.0, # Use all targets (or subset %)
verbose = TRUE
)
## Starting AUC computation...
## Sorting by: expression
## Targets: 100% of available
## Labels: 3 | Cells: 2250
## Backend requested: multisession | backend used: multisession |
## workers: 2
## Note: multisession (SOCK) has worker
## startup/serialization overhead;
## best for larger datasets or Windows compatibility.
##
|
| | 0%
## Loading required namespace: Matrix
##
|
|============ | 25%
## Loading required namespace: Matrix
##
|
|========================= | 50%
## Warning in value[[3L]](cond): Parallel processing failed: BiocParallel errors
## 2 remote errors, element index: 1, 2
## 2 unevaluated and other errors
## first remote error:
## Error in .cell_to_label[[cell_id]]: attempt to select less than one element in get1index
## .
## Falling back to sequential processing...
## AUC computation completed in 14.01 seconds!
## Result: 2250 cells x 10 TFs
## Tip: if runtime is dominated by setup overhead,
## try fewer workers (e.g., 2-4)
## or use backend='multicore' on Linux/macOS.# Extract results in wide format (default)
auc_wide <- get_auc(AS_processor, format = "wide")
head(auc_wide)
## MEF2D E2F4
## AGGGTGATCTGAGGGA-1-NBM-10.138P 0.3451826 0.3714249
## CCTAGCTTCTCCAACC-1-NBM-1.138P 0.3517351 0.3527805
## ATTACTCTCGTGGTCG-1-NBM-10.138P 0.3335810 0.3374751
## TAAGAGAAGCCGCCTA-1-NBM-8.138P 0.3485226 0.3436846
## GATCGATTCAGAGGTG-1-NBM-8.138P 0.3598384 0.3427907
## ACATCAGGTCGCGGTT-1-NBM-10.138P 0.3425284 0.3727721
## SATB1 IRF1
## AGGGTGATCTGAGGGA-1-NBM-10.138P 0.29333333 0.4651859
## CCTAGCTTCTCCAACC-1-NBM-1.138P 0.28000000 0.4530905
## ATTACTCTCGTGGTCG-1-NBM-10.138P 0.21333333 0.4556063
## TAAGAGAAGCCGCCTA-1-NBM-8.138P 0.08666667 0.4462237
## GATCGATTCAGAGGTG-1-NBM-8.138P 0.22000000 0.4508094
## ACATCAGGTCGCGGTT-1-NBM-10.138P 0.20666667 0.4558394
## ATF6 JUN
## AGGGTGATCTGAGGGA-1-NBM-10.138P 0.3616411 0.3961992
## CCTAGCTTCTCCAACC-1-NBM-1.138P 0.3681666 0.4263855
## ATTACTCTCGTGGTCG-1-NBM-10.138P 0.3685420 0.3687974
## TAAGAGAAGCCGCCTA-1-NBM-8.138P 0.3582168 0.4069619
## GATCGATTCAGAGGTG-1-NBM-8.138P 0.3803534 0.3803163
## ACATCAGGTCGCGGTT-1-NBM-10.138P 0.3663562 0.3772032
## JUND POU2F2
## AGGGTGATCTGAGGGA-1-NBM-10.138P 0.3746114 0.3045187
## CCTAGCTTCTCCAACC-1-NBM-1.138P 0.3999677 0.3147107
## ATTACTCTCGTGGTCG-1-NBM-10.138P 0.3817687 0.3028637
## TAAGAGAAGCCGCCTA-1-NBM-8.138P 0.3735578 0.3287177
## GATCGATTCAGAGGTG-1-NBM-8.138P 0.4004624 0.3071023
## ACATCAGGTCGCGGTT-1-NBM-10.138P 0.3786252 0.3140584
## KLF3 BCL11A
## AGGGTGATCTGAGGGA-1-NBM-10.138P 0.2979502 0.2276438
## CCTAGCTTCTCCAACC-1-NBM-1.138P 0.3210761 0.2458956
## ATTACTCTCGTGGTCG-1-NBM-10.138P 0.2975461 0.2211756
## TAAGAGAAGCCGCCTA-1-NBM-8.138P 0.3039077 0.2490154
## GATCGATTCAGAGGTG-1-NBM-8.138P 0.3021116 0.2575864
## ACATCAGGTCGCGGTT-1-NBM-10.138P 0.2894952 0.2296736# Extract results long format
auc_long <- get_auc(AS_processor, format = "long") # Long format
head(auc_long)
## cell tf score label
## 1 AGGGTGATCTGAGGGA-1-NBM-10.138P MEF2D 0.3451826 0
## 2 CCTAGCTTCTCCAACC-1-NBM-1.138P MEF2D 0.3517351 0
## 3 ATTACTCTCGTGGTCG-1-NBM-10.138P MEF2D 0.3335810 0
## 4 TAAGAGAAGCCGCCTA-1-NBM-8.138P MEF2D 0.3485226 0
## 5 GATCGATTCAGAGGTG-1-NBM-8.138P MEF2D 0.3598384 0
## 6 ACATCAGGTCGCGGTT-1-NBM-10.138P MEF2D 0.3425284 0Alternatively you can filter in advance your weights and use SimiCvizExperiment and the wrapper function calculate_activity_scores
weights_df <- read.csv("path/to/your/weights.csv")
weights_df_filtered <- weights_df[weights_df$p_value < 0.01,]
# From CSV files
viz_obj <- SimiCvizExperiment(
weights = weights_df_filtered,
auc = NULL, # Will compute this inthe next section later but can be loaded as well
cell_labels = cell_labels,
label_names = c("control","PD-L1","DAC","Combination"),
colors = c("#e0e0e0", "#a8c8ff", "#ffb6b6", "#c1a9e0"),
meta=list() # Anything you want to store in a list format
)
viz_obj
viz_obj <- calculate_activity_scores(
viz_obj,
expression = expression_mat_path,
adj_r2_threshold = 0.7, # For SimiC style outputs
sort_by="expression", # Rank targets by expression or weight
select_top_k = NULL, # Use all targets (or limit to top K)
percent_of_target = 1.0,
n_cores = 2,
backend = "multisession",
verbose = TRUE
)4.3 Filtering Options
Different GRN methods provide different quality metrics. Adjust filtering based on your method:
# SimiC: Filter by adjusted R² (goodness of fit)
processor_simic <- AUCProcessor(
weights = simic_weights,
expression = expr_mat,
cell_labels = cell_labels,
adj_r2_list = adjusted_r_squared, # a list length as simic_weights
qc_type = "adj_r2",
qc_threshold = 0.7 # Keep targets with R² ≥ 0.7
)
# SCENIC / Pando: Filter by adjusted p-value
processor_scenic <- AUCProcessor(
weights = weights_df,
expression = expr_mat,
cell_labels = cell_labels,
qc_type = "p_value",
qc_threshold = 0.05, # Keep targets with adj_p_val ≤ 0.05
n_cores = 4,
backend = "multisession"
)
# Compute with the same data, different parameters
processor_scenic <- compute_auc(processor_scenic, sort_by = "weight")Once you have computed the activity scores, you are ready to visualize the results!
5 Network Visualization
For visualization we need a SimiCvizExperiment object loaded with the results.
If you used load_SimiCPipeline or computed the activity scores with calculate_activity_scores you already have it.
# Recall above examples
simic_full # Complete SimiCpipeline output
## An object of class SimiCvizExperiment
## 3 label(s), 10 TF(s), 150 target(s)
## Weights: 3 matrices [0: 10 x 150, 1: 10 x 150, 2: 10 x 150]
## AUC: collected (2250 cells x 10 TFs)
## Cell labels: 2250 cells across 3 label(s) [0, 1, 2]
## Label names: 0 = NBM, 1 = SMM, 2 = MM
## Colors: 0 = #3B7EA1, 1 = #E66101, 2 = #B2182B
## TFs: MEF2D, E2F4, SATB1, IRF1, ATF6, JUN, ...
## Targets: RPLP1, RPL36, RPL39, IGHM, EEF1A1, RPS8, ...
## Meta keys: adjusted_r_squared
viz_obj_simic # SimiCPipeline weights -> `calculate_activity_scores`
## An object of class SimiCvizExperiment
## 3 label(s), 10 TF(s), 150 target(s)
## Weights: 3 matrices [0: 10 x 150, 1: 10 x 150, 2: 10 x 150]
## AUC: collected (2250 cells x 10 TFs)
## Cell labels: 2250 cells across 3 label(s) [0, 1, 2]
## Label names: 0 = NBM, 1 = SMM, 2 = MM
## Colors: 0 = #3B7EA1, 1 = #E66101, 2 = #B2182B
## TFs: MEF2D, E2F4, SATB1, IRF1, ATF6, JUN, ...
## Targets: RPLP1, RPL36, RPL39, IGHM, EEF1A1, RPS8, ...
## Meta keys: adjusted_r_squaredOtherwise you can load your results in SimiCvizExperiment as exemplified above including the activity scores results in the auc slot.
Note: Although activity score matrix is not strictly necessary, you need it to fully explore all the visualization tools of SimiCviz.
If you used other methods to calculate cell-specific activity scores with a different algorithm, you can import those as well in:
- wide format: cells x TF scores matrix
- long format: data.frame with cell, tf, score
- .pickle file (only for SimiC v1)
# Create SimiCvizExperiment
simic <- SimiCvizExperiment(weights = simic_weights,
auc = auc_wide,
cell_labels = cell_labels,
label_names = c("NBM","SMM","MM"),
colors = c("#3B7EA1", "#E66101", "#B2182B"),
meta = list(adjusted_r_squared = adjusted_r_squared))
## `weights` provided as a list.
## Cell labels file contains extra columns
## cell / category / label
## AUC in wide format (cells x TFs).
simic
## An object of class SimiCvizExperiment
## 3 label(s), 10 TF(s), 150 target(s)
## Weights: 3 matrices [0: 10 x 150, 1: 10 x 150, 2: 10 x 150]
## AUC: collected (2250 cells x 10 TFs)
## Cell labels: 2250 cells across 3 label(s) [0, 1, 2]
## Label names: 0 = NBM, 1 = SMM, 2 = MM
## Colors: 0 = #3B7EA1, 1 = #E66101, 2 = #B2182B
## TFs: MEF2D, E2F4, SATB1, IRF1, ATF6, JUN, ...
## Targets: RPLP1, RPL36, RPL39, IGHM, EEF1A1, RPS8, ...
## Meta keys: adjusted_r_squaredAfter you created a SimiCvizExperiment object, we can use the plotting functions to visualize the results.
If you want to save your plots in pdf you can define figure output directory out_dir = plot_dir, change the filename and set the argument save = TRUE
plot_dir <- file.path(getwd(),"SimiCviz_output")
dir.create(plot_dir,recursive = TRUE)5.1 Quality Assessment (SimiC only)
To assess the goodness of SimiC’s fitted model of target expression we can plot the distribution of adjusted R² values across all targets.
# Plot distribution of adjusted R² values across targets
# Extract Adjusted R squared from SimiC outputs
out <- read_pickle(weights_file)
adjusted_r_squared <- out$adjusted_r_squared
plot_r2_distribution(adjusted_r_squared, simic, grid = c(1, 3),
save = FALSE, out_dir = plot_dir)
Here we can see that most of the R² values are over 0.7, indicating that the model explains a significant portion of the variance in the target gene expression.
Lower R² might indicate:
- Noisy data
- Missing regulatory interactions
- Complex regulation not captured by linear model
For these targets with low adjusted R2, we suggest to filter them out before plotting. We select a threshold of 0.7 as it was the same used to calculate the activity scores in SimiCPipeline and above. You can change this threshold as needed, but you should re-run the activity score calculation (or SimiCPipeline) with the same threshold to be consistent.
If you already filtered your weight matrix before creating the SimiCviz object you can skip this part
# Select targets by adjusted R2
unselected_targets <- list()
selected_targets <- list()
lab_keys <- names(simic@label_names)
for (lab in lab_keys){
# Save selected for plotting
selected_targets[[lab]] <- simic@target_ids[which(adjusted_r_squared[[lab]] >= 0.7)]
# Save unselected for reporting
label <- simic@label_names[[lab]]
unselected_targets[[label]] <- simic@target_ids[which(adjusted_r_squared[[lab]] < 0.7)]
}
print("Number of unselected targets per label:")
## [1] "Number of unselected targets per label:"
print(sapply(unselected_targets, length))
## NBM SMM MM
## 64 2 15.2 Weights Visualization
5.2.1 TF Barplots:
We will plot the top 30 targets for the top 4 TFs.
You can change the number of TFs and targets to plot by changing the tf_names and top_n arguments.
plot_tf_weights(
simic,
tf_names = simic@tf_ids[1:4],
top_n = 25,
allowed_targets = selected_targets, # Filter by R² if desired
grid = c(2, 2),
save = FALSE,
out_dir = plot_dir,
filename = "TF_weights_barplot.pdf"
)
## Plotting 4 TF weight barplot(s)...
5.2.2 Target Barplots:
Now we will plot the regulators of each target. This plot shows only the TFs with non-zero weights.
plot_target_weights(
simic,
target_names = simic@target_ids[1:4],
labels = c("NBM", "SMM"),
grid = c(2, 2),
save = FALSE,
out_dir = plot_dir,
filename = "Target_weights_barplot.pdf"
)
## Plotting 4 target weight barplot(s)...
Note: If you want to plot all TFs and Targets leave the tf_names and target_names arguments respectively blank.
*Note:: You can also save the output and acces individual plots
all_tfs_barplots <- plot_tf_weights(
simic,
top_n = 25,
grid = NULL,
allowed_targets = selected_targets)all_tfs_barplots[[1]]
5.2.3 Regulatory Network Heatmap
You can extract a TF regulon with the function get_tf_network and plot it as a heatmap.
network <- get_tf_network(simic, "MEF2D", r2_threshold = 0.7)
print(head(network))
## NBM SMM MM
## RPLP1 -0.7501612 0.000000 0.0000000
## RPL36 0.0000000 0.000000 0.0000000
## RPL39 0.0000000 1.230387 0.7876861
## IGHM NA NA 0.0000000
## EEF1A1 0.0000000 0.000000 -1.0652743
## RPS8 0.0000000 0.000000 0.0000000
plot_tf_network_heatmap(simic, "MEF2D",
save = FALSE,
top_n = 15,
r2_threshold = 0.7,
show_values = TRUE,
cmap = c("purple","white","yellow"))
# Dissimilarity Analysis
To prioritize those TFs of interest, that might be driving the regulatory differences between our groups we perform a dissimilarity analysis by comparing per-cell TF activity-score distributions across phenotype labels using the minmax version of total variation distance for multiple distributions. This metric has values between 0 and 1, with values closer to 0 when the group of distributions is more similar to each other. Consequently, higher scores indicate greater regulatory dissimilarity between conditions — i.e., TFs whose activity differs most across phenotypes.
5.3 Global dissimilarity scores
dis_score <- calculate_dissimilarity(simic)
## Calculating dissimilarity scores (all cells)...
## Top 10 TFs by MinMax dissimilarity:
## MinMax_score
## SATB1 0.9360000
## BCL11A 0.9337778
## POU2F2 0.9191111
## IRF1 0.8471111
## E2F4 0.8271111
## KLF3 0.7417778
## JUN 0.6817778
## MEF2D 0.6608889
## JUND 0.6168889
## ATF6 0.5604444
top_tfs <- rownames(dis_score)5.4 Dissimilarity heatmap
We can visualize the top TFs ranked by dissimilarity score in a heatmap fashion.
plot_dissimilarity_heatmap(simic,
top_n = 5,
cmap = "viridis",
save = FALSE)
Sometimes these transcriptional dynamics are cell-type or cell-cluster specific. In this case, we may want to calculate the dissimilarity scores across the labels within specific cell clusters. For that, we can subset the SimiCvizExperiment object by the desired labels and then calculate the dissimilarity scores and plot the heatmap as shown before.
metadata <- read.csv(system.file("extdata/metadata.csv",
package = "SimiCviz"))
# Build cell groups from metadata (e.g. Seurat clusters, cell types, etc.)
cell_groups <- lapply(unique(metadata$cluster),
function(cluster) {
cell_labels$cell[metadata$cluster == cluster]
})
names(cell_groups) <- unique(metadata$cluster)
dissim_grouped <- calculate_dissimilarity(simic, labels = c(1,2),
cell_groups = cell_groups)
## Calculating dissimilarity for selected cell groups...
## Group: healthy (44 cells)
## Group: 17 (19 cells)
## Group: 20 (10 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '20',
## skipping.
## Group: 12 (71 cells)
## Group: 3 (212 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '3',
## skipping.
## Group: 5 (181 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '5',
## skipping.
## Group: 22 (18 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '22',
## skipping.
## Group: 7 (144 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '7',
## skipping.
## Group: 15 (58 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '15',
## skipping.
## Group: 11 (83 cells)
## Group: 19 (14 cells)
## Group: 24 (10 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '24',
## skipping.
## Group: 0 (245 cells)
## Group: 21 (19 cells)
## Group: 9 (91 cells)
## Group: 2 (181 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '2',
## skipping.
## Group: 13 (59 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '13',
## skipping.
## Group: 18 (23 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '18',
## skipping.
## Group: 16 (18 cells)
## Warning in calculate_dissimilarity(simic, labels = c(1, 2),
## cell_groups = cell_groups): <2 labels with cells for '16',
## skipping.
## Top 10 TFs by mean dissimilarity across groups:
## healthy X17 X20 X12 X3 X5 X22 X7 X15
## SATB1 1.000 1.0000000 0 1.0000000 0 0 0 0 0
## POU2F2 1.000 1.0000000 0 1.0000000 0 0 0 0 0
## E2F4 1.000 0.9285714 0 1.0000000 0 0 0 0 0
## MEF2D 0.875 0.8000000 0 0.4853659 0 0 0 0 0
## JUN 0.875 0.6428571 0 0.5650407 0 0 0 0 0
## KLF3 1.000 0.9285714 0 0.5520325 0 0 0 0 0
## IRF1 0.725 0.6428571 0 0.2227642 0 0 0 0 0
## BCL11A 1.000 0.6571429 0 1.0000000 0 0 0 0 0
## ATF6 0.475 0.9285714 0 0.2813008 0 0 0 0 0
## JUND 0.450 0.8571429 0 0.9268293 0 0 0 0 0
## X11 X19 X24 X0 X21
## SATB1 1.0000000 1.0000000 0 1.0000000 1.0000000
## POU2F2 1.0000000 1.0000000 0 1.0000000 0.9333333
## E2F4 0.9830508 1.0000000 0 0.9959016 0.6000000
## MEF2D 0.9491525 1.0000000 0 1.0000000 1.0000000
## JUN 1.0000000 0.9166667 0 1.0000000 0.8000000
## KLF3 1.0000000 1.0000000 0 0.2336066 1.0000000
## IRF1 0.9830508 1.0000000 0 1.0000000 1.0000000
## BCL11A 0.8644068 1.0000000 0 0.2172131 1.0000000
## ATF6 1.0000000 1.0000000 0 0.9180328 0.6833333
## JUND 0.2507062 1.0000000 0 1.0000000 0.7333333
## X9 X2 X13 X18 X16 mean_score
## SATB1 1.0000000 0 0 0 0 0.4210526
## POU2F2 0.9887640 0 0 0 0 0.4169525
## E2F4 0.5955056 0 0 0 0 0.3738437
## MEF2D 0.8539326 0 0 0 0 0.3664974
## JUN 1.0000000 0 0 0 0 0.3578718
## KLF3 0.9775281 0 0 0 0 0.3521968
## IRF1 1.0000000 0 0 0 0 0.3459827
## BCL11A 0.6404494 0 0 0 0 0.3357480
## ATF6 1.0000000 0 0 0 0 0.3308547
## JUND 0.9887640 0 0 0 0 0.3266724
# For all labels
plot_dissimilarity_heatmap(simic,
cell_groups = cell_groups,
top_n = 8,
labels=c(1,2),
cmap=c("magma"),
save = FALSE)
We can also select the labels across which we want to compute the dissimilarity score with the argument labels and customize the color and sorting options.
# For labels 0,2
plot_dissimilarity_heatmap(simic,
cell_groups = cell_groups,
top_n = 5,
sort_by = "healthy",
cmap=c("red", "white", "blue"),
save = FALSE)
6 Activity Score Distributions
Visualize TF activity distributions across conditions. These plots allow you to visually assess the differences in TF activity distributions across conditions, and can help identify TFs with distinct regulatory patterns that may be driving phenotypic differences.
6.0.1 Density plots
This function allows for customization of the plot aesthetics, including fill, transparency, bandwidth adjustment for density estimation, and rug plots to show individual data points. You can also choose to save the plots directly from the function.
# Plot distributions for top TFs
plot_auc_distributions(
simic,
tf_names = top_tfs[1:4],
fill = TRUE,
alpha = 0.6,
bw_adjust = 1/8,
rug = TRUE,
save = FALSE,
out_dir = plot_dir,
filename = "AUC_distributions.pdf",
grid = c(2, 2)
)
# Plot top 4 TFs density distributions
plot_auc_distributions(simic,
labels = c(0,2),
tf_names = top_tfs[1:2],
fill = FALSE,
bw_adjust = 0.5,
rug = FALSE,
out_dir = plot_dir,
filename="AUC_distributions_notfilled_multipage.pdf",
save = FALSE,
grid = c(1,2))
6.0.2 Cumulative distributions (ECDF)
plot_auc_cumulative(
simic,
tf_names = top_tfs[1:4],
rug = TRUE,
grid = c(2, 2),
include_table = TRUE,
save = FALSE,
out_dir = plot_dir
)
6.1 ECDF-based Metrics
Compute area-under-curve metrics from ECDF comparisons:
ecdf_metrics <- calculate_ecdf_auc(simic, tf_names = simic@tf_ids[1:4])
head(ecdf_metrics)
## NBM_ecdf_auc NBM_auc50 NBM_x_at_p50 SMM_ecdf_auc
## MEF2D 0.6486662 0.004245204 0.3484771 0.5888340
## E2F4 0.6407536 0.004853114 0.3598104 0.5649960
## SATB1 0.8176889 0.019600000 0.1733333 0.7013625
## IRF1 0.5447134 0.002492210 0.4554773 0.6456201
## SMM_auc50 SMM_x_at_p50 MM_ecdf_auc MM_auc50
## MEF2D 0.011006497 0.4120845 0.6476606 0.004710800
## E2F4 0.008881756 0.4386415 0.5806568 0.007959580
## SATB1 0.008160962 0.2948023 0.6206071 0.004866559
## IRF1 0.009510367 0.3477034 0.6061756 0.013408928
## MM_x_at_p50 delta_ecdf_auc delta_auc50 delta_x_at_p50
## MEF2D 0.3393122 0.05983225 0.006761294 0.07277230
## E2F4 0.3892008 0.07575763 0.004028642 0.07883117
## SATB1 0.3776027 0.19708184 0.014733441 0.20426940
## IRF1 0.4008829 0.10090673 0.010916718 0.107773926.2 Summary Statistics
6.2.1 Mean activity per TF × phenotype
plot_auc_heatmap(simic, top_n = 20)
6.2.2 Box plots and violin plots
summary_plot <- plot_auc_summary_statistics(simic)
7 Session Info
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets
## [6] methods base
##
## other attached packages:
## [1] SimiCviz_0.99.2 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4
## [3] tidyr_1.3.2 stringi_1.8.7
## [5] lattice_0.22-9 digest_0.6.39
## [7] magrittr_2.0.5 evaluate_1.0.5
## [9] grid_4.6.1 RColorBrewer_1.1-3
## [11] bookdown_0.47 fastmap_1.2.0
## [13] plyr_1.8.9 jsonlite_2.0.0
## [15] Matrix_1.7-6 tinytex_0.60
## [17] gridExtra_2.3.1 BiocManager_1.30.27
## [19] purrr_1.2.2 viridisLite_0.4.3
## [21] scales_1.4.0 codetools_0.2-20
## [23] jquerylib_0.1.4 cli_3.6.6
## [25] rlang_1.3.0 withr_3.0.3
## [27] cachem_1.1.0 yaml_2.3.12
## [29] otel_0.2.0 tools_4.6.1
## [31] parallel_4.6.1 reshape2_1.4.5
## [33] BiocParallel_1.47.0 dplyr_1.2.1
## [35] colorspace_2.1-3 ggplot2_4.0.3
## [37] reticulate_1.46.0 vctrs_0.7.3
## [39] R6_2.6.1 png_0.1-9
## [41] magick_2.9.1 lifecycle_1.0.5
## [43] stringr_1.6.0 pkgconfig_2.0.3
## [45] pillar_1.11.1 bslib_0.11.0
## [47] gtable_0.3.6 glue_1.8.1
## [49] Rcpp_1.1.2 xfun_0.60
## [51] tibble_3.3.1 tidyselect_1.2.1
## [53] knitr_1.51 dichromat_2.0-1
## [55] farver_2.1.2 htmltools_0.5.9
## [57] snow_0.4-4 labeling_0.4.3
## [59] rmarkdown_2.31 compiler_4.6.1
## [61] S7_0.2.2