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.

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.

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")

Loading GRN Data

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.

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 name
target: Target gene name
weight: Edge weight/coefficient
label (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)

SimiCPipeline Outputs

SimiCPipeline is a regularized regression method for inferring condition-specific GRNs. It generates:

.pickle files with TF-target weights per phenotype/label and adjusted R² values for goodness of fit of the target expression model.
.csv file with a TF activity score matrix (cells × TFs)
.pickle files 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.csv

Here 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_squared

If you run the full SimiCPipeline tutorial you can skip Part 2 and go straight to Part 3 for visualizations.

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.

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)
## Downloading uv...Done!
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.000000

Weights from other methods

Most single-cell GRN inference methods produce table outputs with TF-target weights/estimates and p-values.

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 them

Expected 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_val

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     0

The 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))

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  160

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_squared

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.

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 2.91 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.3393418

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 in weights to use (pval, padj, etc.). If ‘qc_type’ == “adj_r2” then cells above qc_threshold will be kept (e.g. R² > 0.7). If ‘qc_type’ != “adj_r2” (we assume is a pvalue) then cells below qc_threshold will 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 6.87 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     0

Alternatively 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
            )

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!

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_squared

Otherwise 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_squared

After 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)

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   1

Weights Visualization

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)...

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]]

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.

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)

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)

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.

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))

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
)

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.10777392

Summary Statistics

Mean activity per TF × phenotype

plot_auc_heatmap(simic, top_n = 20)

Box plots and violin plots

summary_plot <- plot_auc_summary_statistics(simic)

Session Info

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [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: Etc/UTC
## 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] rappdirs_0.3.4      sass_0.4.10        
##  [3] generics_0.1.4      tidyr_1.3.2        
##  [5] stringi_1.8.9       lattice_0.22-9     
##  [7] digest_0.6.39       magrittr_2.0.5     
##  [9] evaluate_1.0.5      grid_4.6.1         
## [11] RColorBrewer_1.1-3  fastmap_1.2.0      
## [13] plyr_1.8.9          rprojroot_2.1.1    
## [15] jsonlite_2.0.0      Matrix_1.7-6       
## [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] here_1.0.2          reticulate_1.46.0  
## [39] buildtools_1.0.0    vctrs_0.7.3        
## [41] R6_2.6.1            png_0.1-9          
## [43] lifecycle_1.0.5     stringr_1.6.0      
## [45] pkgconfig_2.0.3     pillar_1.11.1      
## [47] bslib_0.12.0        gtable_0.3.6       
## [49] glue_1.8.1          Rcpp_1.1.2         
## [51] xfun_0.60           tibble_3.3.1       
## [53] tidyselect_1.2.1    sys_3.4.3          
## [55] knitr_1.51          farver_2.1.2       
## [57] snow_0.4-4          htmltools_0.5.9    
## [59] labeling_0.4.3      rmarkdown_2.31     
## [61] maketools_1.3.2     compiler_4.6.1     
## [63] S7_0.2.2

SimiCviz: Visualization and Analysis of Gene Regulatory Networks