SimiCPipeline - Full Tutorial

Author: Irene Marín-Goñi, PhD student - ML4BM group (CIMA University of Navarra)

This notebook provides a comprehensive guide to how to use the new SimiC pipeline with all available features.

Overview

This tutorial covers: 1. Package installation and setup 2. Complete pipeline initialization with all parameters 3. Running the full pipeline with filtering and AUC calculation steps 4. Follow up steps: Detailed results analysis 5. Analyze previous runs 6. Cross-validation and parameter sweeps

For a simpler introduction, see Tutorial_SimiCPipeline_simple.ipynb

Introduction

SimiC is a GRN inference algorithm for scRNA-Seq data that takes as input: - single-cell imputed expression data (not sparse) - a list of driver genes (transcription factors) - the cell labels (cell phenotypes in order)

and produces TF-specific GRNs for each of the different phenotypes. Given the phenotype order, SimiC adds a similarity constraint when jointly inferring the GRNs for each phenotype, ensuring a smooth transition across phenotypes.

                    0 -> 1; 1 -> 2; 2 -> 3

For more information check our original publication:

Peng, J., Serrano, G., Traniello, I.M. et al. SimiC enables the inference of complex gene regulatory dynamics across cell phenotypes. Commun Biol 5, 351 (2022) DOI: 10.1038/s42003-022-03319-7.

Setup

The easiest way to configure your environment is to follow the README instructions using poetry (or Docker).

Required packages for this tutorial: - simicpipeline - anndata - pandas - numpy - os - pickle

Internally simicpipeline also uses: - scipy - sklearn

We recommend preprocessing the data and set up the directory structure following the Tutorial_SimiCPipeline_preprocessing.ipynb

Import Modules

First, import the necessary modules and set up the path.

import os
print(os.getcwd())
print(os.listdir())
import simicpipeline 
print(f"SimiC pipeline version: ", {simicpipeline.__version__})
from simicpipeline import SimiCPipeline
/home/workdir
['SimiCExampleRun', 'data']
SimiC pipeline version:  {'0.1.0'}

Pipeline steps

Step 1: Initialize the Pipeline

Create a pipeline instance by specifying: - project_dir: Working directory path where input files are located and output files will be saved - run_name: Unique identifier for this analysis run (used as prefix for output files)

We will start this tutorial with the output files generated in the preprocessing tutorial.

print("Initializing SimiC pipeline")
pipeline = SimiCPipeline(
    project_dir="./SimiCExampleRun/KPB25L/Tumor",
    run_name="experiment_tumor"
)
print("\n" + "="*70)
print(f"Pipeline initialized with workdir: {pipeline.project_dir}")
print("="*70)
pipeline.print_project_info()
Initializing SimiC pipeline

======================================================================
Pipeline initialized with workdir: SimiCExampleRun/KPB25L/Tumor
======================================================================
Tumor/
├── inputFiles/
│   ├── TF_list.csv
│   ├── expression_matrix.pickle
│   └── treatment_annotation.csv
└── outputSimic/
    ├── figures/
    └── matrices/
        └── experiment_tumor/

Step 2: Set Input File Paths

Point the pipeline to your input files: - p2df: Path to expression matrix file (genes × cells) stored as a pandas DataFrame in pickle format - p2tf: Path to transcription factor list file (pickle format) containing TF gene names to use as drivers - p2assignment: Path to cell cluster assignment file containing at least one column named label with ordered phenotype labels (integers) matching expression matrix cell order. - df_with_label: Whether the dataframe contains a label column with assignment labels. (Defaut: False)

print("Setting input file paths")
pipeline.set_input_paths(
    p2df = pipeline.project_dir / "inputFiles/expression_matrix.pickle",
    p2tf = pipeline.project_dir / "inputFiles/TF_list.csv",
    p2assignment = pipeline.project_dir / "inputFiles/treatment_annotation.csv"
)
Setting input file paths

Step 3: Set Parameters

Set all available parameters for the SimiC regression:

Core Regularization Parameters:

  • lambda1: L1 regularization strength controlling sparsity in the learned networks (higher values $$ sparser networks with fewer edges, default: \(1e^{-1}\))
  • lambda2: L2 regularization strength controlling similarity between adjacent phenotype networks (higher values $$ more similar networks across phenotypes, default: \(1e^{-5}\))
  • similarity: Boolean flag to enable similarity-based clustering when inferring GRNs (default: True)

Algorithm Parameters:

  • max_rcd_iter: Maximum number of iterations for the coordinate descent optimization algorithm (default: 50000)

IMPORTANT NOTE: This parameter controls the number of iterations of the algorithm and will strongly influence the time it takes to run

Cross-Validation Parameters:

  • cross_val: Boolean flag to enable cross-validation for automatic parameter selection (default: False)
  • k_cross_val: Integer for the number of folds for cross-validation (default: 5)
  • max_rcd_iter_cv: SImilar to max_rcd_iter. To speed up the process the default value is 10000.
  • list_of_l1: List of lambda1 values to test during cross-validation (e.g., [10, 1, 1e-1, 1e-2, 1e-4])
  • list_of_l2: List of lambda2 values to test during cross-validation (e.g., [1e-1, 1e-2, 1e-3, 1e-4, 1e-5])

IMPORTANT NOTE: Cross-validation with multiple lambda values is computationally intensive and may take several hours depending on data size, max_rcd_iter_cv and parameter grid complexity.

print("Setting custom parameters")
pipeline.set_parameters(
    lambda1=1e-1, # Will be optimized if cross-validation = True
    lambda2=1e-2, # Will be optimized by cross-validation = True
    similarity=True,
    max_rcd_iter=5000, # For demonstration purposes we set a smaller number of max iterations, udually 50000
    cross_val=True,
    k_cross_val=4,
    max_rcd_iter_cv=1000, # For demonstration purposes we set a smaller number of max iterations for cross-validation, usually 5000
    list_of_l1=[1e-1, 1e-2],  # Example: [10, 1, 1e-1, 1e-2, 1e-4]
    list_of_l2=[1e-2, 1e-3] # Example: [1e-1, 1e-2, 1e-3, 1e-4, 1e-5]
)
Setting custom parameters

### Step 4: Configure AUC Calculation Parameters

Except for adj_r2_threshold the other parameters are rarely used but are available in SimiC-Jianhao

Define parameters for Area Under Curve calculation which measures TF activity scores: - adj_r2_threshold: Minimum adjusted R² value for filtering target genes based on regression quality (range: 0-1, typical: 0.7) - select_top_k_targets: Number of top-ranked targets to include per TF (None: use all targets, integer: specific number) - percent_of_target: Percentage of targets to include after ranking (range: 0-1, default: 1 for 100%) - sort_by: Criterion for ranking and selecting targets (‘expression’: by mean expression level, ‘weight’: by network weight, ‘adj_r2’: by adjusted R² value) - num_cores: Number of CPU cores to use for parallel processing (-1: use all available cores, numbers below -1 will use [#total cores + 1 + num_cores]). Default: 1.

auc_params = {
    'adj_r2_threshold': 0.7,
    'select_top_k_targets': None,
    'percent_of_target': 1,
    'sort_by': 'expression',
    'num_cores': -2
}

Validate the inputs before run the pipeline

pipeline.validate_inputs()
✓ All required input files found

Step 5: Run the Complete Pipeline

Executing the full pipeline in one command is easy with this implementation.

Steps incuded: 1. Input validation 2. SimiC regression 3. Weight filtering 4. AUC calculation

Custom Pipeline Execution Options:

  • skip_filtering: If True, skips the weight filtering step that removes low-importance edges (default: False)
  • calculate_raw_auc: If True, calculates AUC scores on unfiltered regression weights before any filtering (default: False)
  • calculate_filtered_auc: If True, calculates AUC scores on filtered weights after removing noise (default: True)
  • variance_threshold: Threshold for TF selection based on variance explained - keeps TFs explaining this cumulative percentage of target variance (range: 0-1, common value: 0.9 for 90%)
  • auc_params: Dictionary containing all AUC calculation parameters as defined in Step4
print("Running complete SimiC pipeline...")
print("This may take several minutes depending on data size.\n")

pipeline.run_pipeline(
    skip_filtering=False,
    calculate_raw_auc=False,
    calculate_filtered_auc=True,
    variance_threshold=0.9,
    auc_params=auc_params
)
Running complete SimiC pipeline...
This may take several minutes depending on data size.


======================================================================
SIMIC PIPELINE
======================================================================

==================================================
Running SimiC Regression
Run name: experiment_tumor
✓ All required input files found
Running cross-validation with following lambdas: [0.1, 0.01] (L1), [0.01, 0.001] (L2)
==================================================

Expression matrix for regression shape =  (21490, 1100)
df test =  (4298, 1100)
test data assignment set: {0, 1, 2, 3}
df train =  (17192, 1100)
train data assignment set: {0, 1, 2, 3}
-------
-------
.... generating train set
cell type  0
    TF size: (5034, 101)
    Target size: (5034, 1000)
cell type  1
    TF size: (3929, 101)
    Target size: (3929, 1000)
cell type  2
    TF size: (4754, 101)
    Target size: (4754, 1000)
cell type  3
    TF size: (3475, 101)
    Target size: (3475, 1000)
-------
.... generating test set
cell type  0
    TF size: (1205, 101)
    Target size: (1205, 1000)
cell type  1
    TF size: (1035, 101)
    Target size: (1035, 1000)
cell type  2
    TF size: (1141, 101)
    Target size: (1141, 1000)
cell type  3
    TF size: (917, 101)
    Target size: (917, 1000)
-------

----------------------------------------------------------------------
Start cross validation!!!
----------------------------------------------------------------------

Trying lambda1 = [0.1, 0.01] and lambda2 = [0.01, 0.001]
lambda1 = 0.1, lamda2 = 0.01, done
----> Averaged adjusetd R2 = 0.4244
----> R2 in folds = [0.4512, 0.4380, 0.4167, 0.3918]
////////
lambda1 = 0.1, lamda2 = 0.001, done
----> Averaged adjusetd R2 = 0.4243
----> R2 in folds = [0.4514, 0.4379, 0.4163, 0.3916]
////////
lambda1 = 0.01, lamda2 = 0.01, done
----> Averaged adjusetd R2 = 0.6407
----> R2 in folds = [0.6647, 0.6571, 0.6355, 0.6054]
////////
lambda1 = 0.01, lamda2 = 0.001, done
----> Averaged adjusetd R2 = 0.6411
----> R2 in folds = [0.6651, 0.6555, 0.6374, 0.6065]
////////
Cross Validation done! 

Selected: lambda1 = 0.01, lambda2 = 0.001, opt R squared on eval 0.6411
----------------------------------------------------------------------

----------------------------------------------------------------------
Begin the SimiC Regression!!!
----------------------------------------------------------------------

-------
final train error w.o. reg = 203.9880 +/- 0.0000 SD
test error w.o. reg = 210.6154 +/- 0.0000 SD
-------
R squared of test set (before): -14.5708
R squared of train set (after): 0.8148 +/- 0.0000 SD
R squared of test set (after): 0.8039 +/- 0.0000 SD
Pickle method succeeded, saved to SimiCExampleRun/KPB25L/Tumor/outputSimic/matrices/experiment_tumor/experiment_tumor_L1_0.01_L2_0.001_simic_matrices.pickle
Updated lambdas from cross-validation: L1=0.01, L2=0.001

✓ SimiC regression completed in 5min 37s

==================================================
Filtering weights using BIC criterion
Variance threshold: 0.9
==================================================

Loaded weights for 4 phenotype labels
Number of TFs: 100
Number of targets: 1000

Processing label 0...
  Targets with non-zero weights: 1000/1000
  TFs kept per target: Mean=21.86, Median=22, Max=40, Min=7

Processing label 1...
  Targets with non-zero weights: 1000/1000
  TFs kept per target: Mean=23.06, Median=24, Max=43, Min=8

Processing label 2...
  Targets with non-zero weights: 1000/1000
  TFs kept per target: Mean=22.32, Median=23, Max=41, Min=7

Processing label 3...
  Targets with non-zero weights: 1000/1000
  TFs kept per target: Mean=21.00, Median=21, Max=41, Min=2

✓ Weight filtering completed in 32s
Filtered weights saved to: SimiCExampleRun/KPB25L/Tumor/outputSimic/matrices/experiment_tumor/experiment_tumor_L1_0.01_L2_0.001_simic_matrices_filtered_BIC.pickle

==================================================
Calculating AUC matrices (filtered weights)
==================================================

Using all available cores for parallel processing.
Pickle method succeeded, saved to SimiCExampleRun/KPB25L/Tumor/outputSimic/matrices/experiment_tumor/experiment_tumor_L1_0.01_L2_0.001_wAUC_matrices_filtered_BIC.pickle

✓ AUC calculation completed in 9min 33s

Collecting AUC for all labels...
✓ Collected AUC for all labels saved to: SimiCExampleRun/KPB25L/Tumor/outputSimic/matrices/experiment_tumor/experiment_tumor_L1_0.01_L2_0.001_wAUC_matrices_filtered_BIC_collected.csv

======================================================================
PIPELINE EXECUTION SUMMARY
======================================================================

Run name: experiment_tumor
Project directory: SimiCExampleRun/KPB25L/Tumor

Parameters:
  - Lambda1: 0.01
  - Lambda2: 0.001
  - Number of TFs: 100
  - Number of targets: 1000

Timing:
  - simic_regression: 5min 37s
  - filtering: 32s
  - auc_filtered: 9min 33s
  - total: 15min 46s

Available Results:
✓ Ws_raw
✓ Ws_filtered
✗ auc_raw
✓ auc_filtered

======================================================================

======================================================================

Success!

pipeline.print_project_info(max_depth=4)
Tumor/
├── inputFiles/
│   ├── TF_list.csv
│   ├── expression_matrix.pickle
│   └── treatment_annotation.csv
└── outputSimic/
    ├── figures/
    └── matrices/
        └── experiment_tumor/
            ├── experiment_tumor_L1_0.01_L2_0.001_simic_matrices.pickle
            ├── experiment_tumor_L1_0.01_L2_0.001_simic_matrices_filtered_BIC.pickle
            ├── experiment_tumor_L1_0.01_L2_0.001_wAUC_matrices_filtered_BIC.pickle
            └── experiment_tumor_L1_0.01_L2_0.001_wAUC_matrices_filtered_BIC_collected.csv

How to continue?

In this section we will show a quick overview of all results generated by the SimiCPipeline run. All the important results are automatically saved but if you want to explore and save some of them in different format here are some examples on how to do that.

For advanced visualizations see the Tutorial_SimiCPipeline_visualization

Check Available Results

If you stoped the pipeline and want to resume it later, you should re-initialize the pipeline using the same paths and lambda parameters. This way pipeline.available_results() will look for filenames in the corresponding directory and check if they exist.

import pandas as pd
import numpy as np
from simicpipeline import SimiCPipeline

pipeline = SimiCPipeline(
    project_dir="./SimiCExampleRun/KPB25L/Tumor",
    run_name="experiment_tumor"
)
pipeline.set_input_paths(
    p2df = pipeline.project_dir / "inputFiles/expression_matrix.pickle",
    p2assignment = pipeline.project_dir / "inputFiles/treatment_annotation.csv",
    p2tf = pipeline.project_dir / "inputFiles/TF_list.csv"
)
pipeline.set_parameters(
    lambda1=1e-2,
    lambda2=1e-3
)
pipeline._print_summary()

======================================================================
PIPELINE EXECUTION SUMMARY
======================================================================

Run name: experiment_tumor
Project directory: SimiCExampleRun/KPB25L/Tumor

Parameters:
  - Lambda1: 0.01
  - Lambda2: 0.001
  - Number of TFs: 100
  - Number of targets: 1000

Timing:

Available Results:
✓ Ws_raw
✓ Ws_filtered
✗ auc_raw
✓ auc_filtered

======================================================================

======================================================================

Load and Inspect Results

Here we show how to access each item generated in the pipeline for your convinience (you can further analyze/explore them).

We designed a load_results function so you can easily access all the pipeline results. - result_type: can be 'Ws_raw', 'Ws_filtered', 'auc_raw' or 'auc_filtered'

Load Filtered Weights

Access the filtered weight matrices for downstream analysis.

print("="*70)
print("Accessing filtered weights...\n")

simic_results = pipeline.load_results('Ws_filtered')
print(f"SimiC results keys: {list(simic_results.keys())}")

weight_dic = simic_results['weight_dic']
print(f"\nWeight dictionary keys (labels): {list(weight_dic.keys())}")
print(f"Number of cell populations: {len(weight_dic)}")

# Inspect first population
first_label = list(weight_dic.keys())[0]
first_weights = weight_dic[first_label]
print(f"\nWeight matrix for label {first_label}:")
print(f"  Shape: {first_weights.shape}")
print(f"  Non-zero entries: {(first_weights != 0).sum()}")
pd.DataFrame(first_weights[0:3,0:5],index = simic_results['TF_ids'][0:3], columns=simic_results['query_targets'][0:5])
======================================================================
Accessing filtered weights...

SimiC results keys: ['weight_dic', 'adjusted_r_squared', 'standard_error', 'TF_ids', 'query_targets']

Weight dictionary keys (labels): [0, 1, 2, 3]
Number of cell populations: 4

Weight matrix for label 0:
  Shape: (101, 1000)
  Non-zero entries: 22858
Malat1 Rn18s-rs5 Cp Brinp3 Cmss1
Pms1 0.0 0.0 0.000000 0.0 0.0
Tshz1 0.0 0.0 -0.919577 0.0 0.0
Wdhd1 0.0 0.0 0.000000 0.0 0.0

Analyze Weight Distribution

Generate summary statistics of the learned weights.

This function compares the sparsity of the network before and after filtering weights by BIC criterion (keeping only the top TFs that explain at least variance_threshold of the variance for each target gene, Step5).

pipeline.analyze_weights()

======================================================================
ANALYZING WEIGHT MATRICES
======================================================================

Label 0:
  Raw weights:
    - Sparsity: 0.00%
    - Avg non-zero TFs per target: 100.00
  Filtered weights:
    - Sparsity: 78.14%
    - Avg non-zero TFs per target: 21.86
  Reduction: 78.14% more sparse

Label 1:
  Raw weights:
    - Sparsity: 0.00%
    - Avg non-zero TFs per target: 100.00
  Filtered weights:
    - Sparsity: 76.94%
    - Avg non-zero TFs per target: 23.06
  Reduction: 76.94% more sparse

Label 2:
  Raw weights:
    - Sparsity: 0.00%
    - Avg non-zero TFs per target: 100.00
  Filtered weights:
    - Sparsity: 77.68%
    - Avg non-zero TFs per target: 22.32
  Reduction: 77.68% more sparse

Label 3:
  Raw weights:
    - Sparsity: 0.00%
    - Avg non-zero TFs per target: 100.00
  Filtered weights:
    - Sparsity: 79.00%
    - Avg non-zero TFs per target: 21.00
  Reduction: 79.00% more sparse

Load Filtered AUC Scores

Examine the filtered AUC scores.

print("="*70)
print("Accessing filtered AUC scores...\n")

auc_filtered = pipeline.load_results('auc_filtered')
print(f"AUC dictionary keys (labels): {list(auc_filtered.keys())}")
print(f"Number of cell populations: {len(auc_filtered)}")
print("Shape of AUC matrix for first label (unprocessed):", auc_filtered[0].shape)
print("="*70)
======================================================================
Accessing filtered AUC scores...

AUC dictionary keys (labels): [0, 1, 2, 3]
Number of cell populations: 4
Shape of AUC matrix for first label (unprocessed): (21490, 100)
======================================================================

Note that each key in the dictionary corresponds to a label in the data but the AUC matrices outputted have values for all the cells (original_matrix shape)

You need to subset the AUC matrix to get only the cells corresponding to that label.

For that we have the following handy function (also used internally by pipeline.subset_label_specific_auc() that loads the AUC data and subsets for specific cells as shown below:

auc_subset = pipeline.subset_label_specific_auc( result_type = 'auc_filtered',label=0)
# you can now calculate the statistics on the subsetted AUC matrix
import numpy as np
print(f"  AUC score statistics:")
print(f"  Shape: {auc_subset.shape} (cells x TFs)")
print(f"    - Mean: {np.nanmean(auc_subset.values):.4f}")
auc_subset.head(3)
  AUC score statistics:
  Shape: (6239, 100) (cells x TFs)
    - Mean: 0.4592
Pms1 Tshz1 Wdhd1 Rere Tead4 Mecom Nr1h5 Tead1 Foxn3 Tfdp2 ... Creb3l2 Zfp950 Trps1 Hivep2 Aff2 Gtf2i Bbx Fosl1 Arid5b Tulp4
01_01_28__s1 0.401697 0.331343 0.312967 0.457909 0.419254 0.335233 0.398261 0.484048 0.439490 0.452826 ... 0.255514 0.518625 0.438519 0.539573 0.418137 0.397457 0.611732 0.640829 0.439133 0.518133
01_01_62__s1 0.301913 0.361450 0.211921 0.532261 0.418367 0.422623 0.285906 0.468733 0.481279 0.393999 ... 0.322578 0.511813 0.531098 0.665322 0.494860 0.594714 0.725136 0.462431 0.611774 0.663357
01_02_38__s1 0.307675 0.389989 0.227622 0.521038 0.448241 0.429574 0.280809 0.434074 0.541557 0.462338 ... 0.430952 0.503680 0.489732 0.625177 0.485962 0.504431 0.573783 0.537172 0.502197 0.556878

3 rows × 100 columns

Notice in the output of Step 5 that we have generated a file named <id_name>_wAUC_matrices_filtered_BIC_collected.csv. This data frame contains the activity scores already extracted for each of your phenotypes. You can easily incorporate this data into your Seurat/SingleCellExperiment/AnnData metadata and visualize it in UMAP. See Tutorial_SImiCPipeline_visualization.ipynb for an examples.

If you want to access all the activity scores you can just reload the data frame:

import pandas as pd
# Complete File name SimiCExampleRun/KPB25L/Tumor/outputSimic/matrices/experiment_tumor/experiment_tumor_L1_0.1_L2_0.01_wAUC_matrices_filtered_BIC_collected.csv
auc_subset_all = pd.read_csv(pipeline.p2auc_filtered.with_name(pipeline.p2auc_filtered.stem + "_collected.csv"), index_col=0)
print("Shape of collected AUC subset:", auc_subset_all.shape)
auc_subset_all.head(3)
Shape of collected AUC subset: (21490, 100)
Pms1 Tshz1 Wdhd1 Rere Tead4 Mecom Nr1h5 Tead1 Foxn3 Tfdp2 ... Creb3l2 Zfp950 Trps1 Hivep2 Aff2 Gtf2i Bbx Fosl1 Arid5b Tulp4
01_01_28__s1 0.401697 0.331343 0.312967 0.457909 0.419254 0.335233 0.398261 0.484048 0.439490 0.452826 ... 0.255514 0.518625 0.438519 0.539573 0.418137 0.397457 0.611732 0.640829 0.439133 0.518133
01_01_62__s1 0.301913 0.361450 0.211921 0.532261 0.418367 0.422623 0.285906 0.468733 0.481279 0.393999 ... 0.322578 0.511813 0.531098 0.665322 0.494860 0.594714 0.725136 0.462431 0.611774 0.663357
01_02_38__s1 0.307675 0.389989 0.227622 0.521038 0.448241 0.429574 0.280809 0.434074 0.541557 0.462338 ... 0.430952 0.503680 0.489732 0.625177 0.485962 0.504431 0.573783 0.537172 0.502197 0.556878

3 rows × 100 columns

Or re-generated it like this:

import pandas as pd
auc_subset_list = []
for label in [0,1,2,3]:
    auc_subset = pipeline.subset_label_specific_auc('auc_filtered', label=label)
    auc_subset_list.append(auc_subset)
auc_subset_all = pd.concat(auc_subset_list, axis=0)
print("Shape of collected AUC subset:", auc_subset_all.shape)
auc_subset_all.head(3)
# auc_subset_all.to_csv("auc_filtered_collected.csv")
Shape of collected AUC subset: (21490, 100)
Pms1 Tshz1 Wdhd1 Rere Tead4 Mecom Nr1h5 Tead1 Foxn3 Tfdp2 ... Creb3l2 Zfp950 Trps1 Hivep2 Aff2 Gtf2i Bbx Fosl1 Arid5b Tulp4
01_01_28__s1 0.401697 0.331343 0.312967 0.457909 0.419254 0.335233 0.398261 0.484048 0.439490 0.452826 ... 0.255514 0.518625 0.438519 0.539573 0.418137 0.397457 0.611732 0.640829 0.439133 0.518133
01_01_62__s1 0.301913 0.361450 0.211921 0.532261 0.418367 0.422623 0.285906 0.468733 0.481279 0.393999 ... 0.322578 0.511813 0.531098 0.665322 0.494860 0.594714 0.725136 0.462431 0.611774 0.663357
01_02_38__s1 0.307675 0.389989 0.227622 0.521038 0.448241 0.429574 0.280809 0.434074 0.541557 0.462338 ... 0.430952 0.503680 0.489732 0.625177 0.485962 0.504431 0.573783 0.537172 0.502197 0.556878

3 rows × 100 columns

Analyze AUC Scores

Examine the distribution of TF-target predicted activity scores in each phenotype label.

pipeline.analyze_auc_scores()

======================================================================
ANALYZING AUC SCORES
======================================================================

Label 0:
  Shape: (6239, 100) (cells x TFs)
  AUC score statistics:
    - Mean: 0.4592
    - Median: 0.4707
    - Std: 0.1196
    - Min: 0.1186
    - Max: 0.9230
  Top 5 TFs by average AUC:
    - Bbx: 0.6251
    - Hivep2: 0.6143
    - Arnt: 0.6104
    - Baz2b: 0.6059
    - Nfia: 0.6030

Label 1:
  Shape: (4964, 100) (cells x TFs)
  AUC score statistics:
    - Mean: 0.4519
    - Median: 0.4598
    - Std: 0.1149
    - Min: 0.1120
    - Max: 0.9788
  Top 5 TFs by average AUC:
    - Tcf12: 0.6197
    - Ncor2: 0.6185
    - Rfx3: 0.6145
    - Tulp4: 0.6065
    - Bbx: 0.5929

Label 2:
  Shape: (5895, 100) (cells x TFs)
  AUC score statistics:
    - Mean: 0.4372
    - Median: 0.4443
    - Std: 0.1220
    - Min: 0.0270
    - Max: 0.8669
  Top 5 TFs by average AUC:
    - Baz2b: 0.6936
    - Rere: 0.6284
    - Mef2a: 0.6115
    - Zmiz1: 0.5933
    - Tcf12: 0.5803

Label 3:
  Shape: (4392, 100) (cells x TFs)
  AUC score statistics:
    - Mean: 0.4238
    - Median: 0.4342
    - Std: 0.1155
    - Min: 0.0060
    - Max: 0.9749
  Top 5 TFs by average AUC:
    - Tcf12: 0.6564
    - Rere: 0.6394
    - Runx1: 0.5713
    - Baz2b: 0.5670
    - Cux1: 0.5511

Calculate Dissimilarity Scores

Compute dissimilar regulons between phenotype populations.

All Labels

If all labels are selected it will calculate the min-max difference among them (independently of the label order: 0,1,2)

print("----> For all labels:")
MinMax_all = pipeline.calculate_dissimilarity()
print(f"\nDissimilarity matrix shape: {MinMax_all.shape}")
print(type(MinMax_all))
# MinMax_all.head()
----> For all labels:

======================================================================
CALCULATING DISSIMILARITY SCORES ACROSS LABELS
======================================================================


Calculating dissimilarity scores...

Top 10 TFs by MinMax dissimilarity score:
  Zfp950: 0.6234
  Nfat5: 0.6198
  Twist2: 0.6087
  Esr1: 0.6077
  Tead1: 0.5861
  Etv6: 0.5816
  Runx1: 0.5729
  Tcf7l2: 0.5686
  Arid1b: 0.5686
  Trps1: 0.5588

Dissimilarity matrix shape: (100, 1)
<class 'pandas.core.frame.DataFrame'>

Specific Labels

Calculate dissimilarity between selected labels only. - select_labels: List of integer labels to include in dissimilarity calculation (e.g., [0, 3] for comparing populations 0 and 3)

MinMax_0_3 = pipeline.calculate_dissimilarity(select_labels=[0, 3])

======================================================================
CALCULATING DISSIMILARITY SCORES ACROSS LABELS
======================================================================

Comparing labels [0, 3].

Calculating dissimilarity scores...

Top 10 TFs by MinMax dissimilarity score:
  Runx1: 0.9709
  Tead1: 0.9571
  Esr1: 0.9453
  Rere: 0.9186
  Zfp950: 0.9047
  Etv5: 0.9005
  Trps1: 0.8617
  Ets1: 0.8521
  Nfia: 0.8162
  Etv6: 0.8145
  • result_type: Name of the AUC results to load (‘auc_raw’ for unfiltered or ‘auc_filtered’ for filtered weights)
  • label: Integer specifying which cell phenotype/population to extract (must match labels in assignment file)

Extract TF-Specific Networks

Get the regulatory network for specific transcription factors. - TF_name: Name of the transcription factor gene (must be present in the TF list provided to the pipeline) - stacked: If True, returns a pandas DataFrame with GRN weights for all labels in separate columns; if False, returns a dictionary of separate pandas Series per label

# Example: Get network for Bnc2
tf_name = "Bnc2"
network = pipeline.get_TF_network(TF_name= tf_name, stacked=True)

print(f"Network for {tf_name}:")
print(f"  Shape: {network.shape}")
print(f"  Columns: {list(network.columns)}")
print(f"\nTop 10 targets:")
print(network.head(10))
Retrieving network for TF: Bnc2
Network for Bnc2:
  Shape: (809, 4)
  Columns: [0, 1, 2, 3]

Top 10 targets:
                  0         1         2         3
Malat1     0.000000  1.534094  0.789632  1.437999
Rn18s-rs5 -1.547416  0.000000  0.000000 -1.029025
Cp         0.974170 -1.060997  0.915877  1.299228
Brinp3     2.762555  3.653162  4.099088  3.709227
Cmss1     -0.939836  0.000000  0.000000 -0.881202
Pcdh7      0.000000  1.529311  1.393673  4.440386
Pnrc1      1.777511 -1.833807  1.352006  0.000000
Sema5a     2.843518  2.130922  2.650287  3.768668
Slit2      0.000000  2.223871  5.097978  5.389786
Cdk14      3.651989  2.647673  4.229523  3.920700
# Example: Get auc for Bnc2
tf_name = ["Bnc2"]
tf_auc = pipeline.get_TF_auc(TF_name = tf_name, stacked=False)
print(type(tf_auc),"\n")

tf_auc[0].head()
Retrieving AUC for TF: ['Bnc2']
<class 'dict'>
Bnc2
01_01_28__s1 0.381035
01_01_62__s1 0.473624
01_02_38__s1 0.433434
01_02_54__s1 0.389323
01_02_81__s1 0.477805 
tf_auc_stack = pipeline.get_TF_auc(TF_name= tf_name, stacked=True)
print(type(tf_auc_stack),"\n")
tf_auc_stack.head()
Retrieving AUC for TF: ['Bnc2']
<class 'pandas.core.frame.DataFrame'>
Bnc2 label category
01_01_28__s1 0.381035 0 control
01_01_62__s1 0.473624 0 control
01_02_38__s1 0.433434 0 control
01_02_54__s1 0.389323 0 control
01_02_81__s1 0.477805 0 control 

pipeline.p2simic_matrices
PosixPath('SimiCExampleRun/KPB25L/Tumor/outputSimic/matrices/experiment_tumor/experiment_tumor_L1_0.01_L2_0.001_simic_matrices.pickle')
tf_names = ["Bnc2", "Runx2"]
tf_auc = pipeline.get_TF_auc(TF_name= tf_names, stacked=True)
tf_auc.head()
Retrieving AUC for TF: ['Bnc2', 'Runx2']
Bnc2 Runx2 label category
01_01_28__s1 0.381035 0.365570 0 control
01_01_62__s1 0.473624 0.531623 0 control
01_02_38__s1 0.433434 0.489543 0 control
01_02_54__s1 0.389323 0.377849 0 control
01_02_81__s1 0.477805 0.516053 0 control 
Label 2 processed in 557.52 seconds
Label 1 processed in 558.30 seconds
Label 3 processed in 566.28 seconds
Label 0 processed in 570.77 seconds

Summary

This tutorial covered:

✓ Complete pipeline initialization with all parameters

✓ Running the full analysis with filtering and AUC calculation

✓ Loading and inspecting all result types

✓ Analyzing weights and AUC score distributions

✓ Calculating dissimilarity between cell populations

✓ Extracting TF-specific regulatory networks

Next Steps

  1. Visualize the results: Check Tutorial_SimiCPipeline_visualization.ipynb
  2. Optimize parameters: Cross-validation help you to find optimal λ₁ and λ₂ values for your data but sometimes you need to re-assess.
  3. Explore results: Investigate specific TFs of interest and their predicted target networks
  4. Validate findings: Compare inferred regulatory relationships with known literature or experimental data
  5. Visualize networks: Create network diagrams for key TFs using network visualization tools

Additional Resources

  • Data preprocessing tutorial: Tutorial_SimiCPipeline_preprocessing.ipynb
  • Visualization tutorial: Tutorial_SimiCPipeline_visualization.ipynb
  • Check example (non interactive) scripts: run_preprocessing.py and run_simicpipeline.py
  • API documentation: Check the SimiCPipeline class documentation for detailed information
  • Original publication: Peng et al., Commun Biol 5, 351 (2022)
  • High impact publications using SimiC algorithm: