SimiC Preprocessing Pipeline Tutorial

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

This notebook demonstrates how to preprocess single-cell RNA-seq data for SimiC analysis.

Overview

This preprocessing tutorial covers: 1. Package installation and setup 2. MAGIC imputation pipeline 3. Gene selection and experiment setup 4. Preparing input files for SimiC

For running SimiC analysis see Tutorial_SimiCPipeline_full.ipynb

Introduction

Before running SimiC, you need to: - Impute your scRNA-seq data. We recommend to use MAGIC and include a wrapper class MagicPipeline to ease the process. - Select top variable genes based on Median Absolute Deviation (MAD) or the genes of interest from which you want to infer the gene regulatory network. - Prepare input files in the correct format for SimiCPipeline.

This tutorial shows you how to do all of this using the SimiC preprocessing modules.

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 - scprep (in preprocessing) - magic-impute (in preprocessing)

Import Modules

First, import the necessary preprocessing modules.

import os
os.chdir('/home/workdir/')
print(os.getcwd())
print(os.listdir())
import simicpipeline 
print(f"SimiC pipeline version: ", {simicpipeline.__version__})

/home/workdir
['cell_ids.txt', 'data', 'SIMIC_bv', 'run_simicpipeline-CS-NBM-MM.log', 'CaseStudy', 'wnn_noMM2_processed.rds', 'metadata.csv', 'NBM-MM', 'multiome-final-data', 'Boiarsky', '.ipynb_checkpoints', 'Rplots.pdf', 'data_GSE145977', 'Dockerfile', 'genes_ids.txt', 'singlecell_matrix.mtx', 'SIMIC_RUN', 'cell_ids_bv.txt', 'SIMIC_RUN2', 'run_simicpipeline-CS-NBM-SMM-MM.log', '.claude', 'run_simicpipeline-CS.log', 'SimicViz', '.Trash-0', 'test1.ipynb', 'metadata_bv.csv', 'SIMIC_RUN3', 'testing.ipynb', 'genes_ids_bv.txt', 'run_simicpipeline-CS-MT.log', 'singlecell_matrix_bv.mtx', 'NBM-SMM-MM']
SimiC pipeline version:  {'0.1.0'}

# Part 1: MAGIC Imputation Pipeline

MAGIC (Markov Affinity-based Graph Imputation of Cells) is used to denoise and impute scRNA-seq data. MagicPipeline facilitates the steps described in Magic Tutorial

Step 1.1: Load Your Data

Load your raw expression data. Note that MacigPipeline class expects AnnData format.

Below you will see different examples on how to generate the AnnData object from different input files (including Seurat if you are more familiar with R)

Note 1: If you have already filtered and transformed your data but want to repeat following these steps, make sure the adata object has the raw counts in the adata.raw.X slot

Note 2: If you have already processed (inputed) your data you can jump to Part 2 Experiment Setup

From Seurat

If you have your data in a Seurat object (R package) the easiest approach is this:

write.table(data.frame(Cells = colnames(seurat_obj)), file.path("/path/to/data/", "cell_ids.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)

write.table(data.frame(Genes = rownames(seurat_obj)), file.path("/path/to/data/", "genes_ids.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)

write.table(seurat_obj@metadata, file.path("/path/to/data/", "metadata.csv"), sep = ",",row.names = TRUE, col.names = TRUE, quote = FALSE)

m_raw = GetAssayData(seurat_obj, assay = "RNA", layer = "counts") # Take the raw counts

Matrix::writeMM(m_raw, paste0(magic_path, "/singlecell_matrix.mtx")) # Save it in MatrixMarket format 

Then follow the example code below.

print("Load your AnnData object here")

# # Example: Load from Matrix Market format
import pandas as pd
import anndata as ad
from pathlib import Path

# This will return a pd.DataFrame
df = simicpipeline.load_from_matrix_market( 
    matrix_path=Path("./Boiarsky/output/singlecell_matrix.mtx"),
    genes_path=Path("./Boiarsky/output/genes_ids.txt"),
    cells_path=Path("./Boiarsky/output/cell_ids.txt"),
    transpose=False, # from python! 
    cells_index_name="Cell",
)

adata = ad.AnnData(X=df.values, obs=pd.DataFrame(index=df.index), var=pd.DataFrame(index=df.columns))

obs_meta = pd.read_csv("./Boiarsky/output/metadata.csv", sep=",", index_col=0)
print(obs_meta.shape)
print(obs_meta.head())

# Match the obs metadata index to the adata.obs_names
obs_meta = obs_meta.loc[adata.obs_names]
adata = ad.AnnData(adata.X, obs=obs_meta, var=adata.var)
# If your data is raw, you should set it properly in AnnData object
adata.raw = adata.copy()

Load your AnnData object here
(28550, 7)
                             sample_ID disease_stage  n_genes  frac_mito  \
index                                                                      
AAACCTGAGGTAAACT-1-MM-5.138P      MM-5            MM     1151   0.129362   
AAACCTGGTCCCTACT-1-MM-5.138P      MM-5            MM     1260   0.108620   
AAACCTGGTTACCGAT-1-MM-5.138P      MM-5            MM     1347   0.058278   
AAACGGGAGCACACAG-1-MM-5.138P      MM-5            MM     1941   0.032493   
AAACGGGTCACAATGC-1-MM-5.138P      MM-5            MM      889   0.071260   

                              n_counts cluster normal_or_neoplastic  
index                                                                
AAACCTGAGGTAAACT-1-MM-5.138P    5875.0       9           neoplastic  
AAACCTGGTCCCTACT-1-MM-5.138P    6021.0       9           neoplastic  
AAACCTGGTTACCGAT-1-MM-5.138P    8082.0       9           neoplastic  
AAACGGGAGCACACAG-1-MM-5.138P   17296.0       9           neoplastic  
AAACGGGTCACAATGC-1-MM-5.138P    3817.0       9           neoplastic  

adata.obs["disease_stage"] = pd.Categorical(
    adata.obs["disease_stage"],
    categories=["NBM", "SMM", "MM"],
    ordered=True
)

# Check it worked
print(adata.obs["disease_stage"].value_counts())
print(adata.obs["disease_stage"].cat.categories)

disease_stage
MM     10790
NBM     9329
SMM     8431
Name: count, dtype: int64
Index(['NBM', 'SMM', 'MM'], dtype='object')

adata

AnnData object with n_obs × n_vars = 28550 × 22273
    obs: 'sample_ID', 'disease_stage', 'n_genes', 'frac_mito', 'n_counts', 'cluster', 'normal_or_neoplastic'

Step 1.2: Initialize MAGIC Pipeline

Create a MAGIC pipeline instance: - input_data: Your AnnData object. If you run the full pipline starting from raw counts they should be in adata.raw.X - project_dir: Project directory where magic_outputdir will be created and files will be saved - magic_output_file: Filename for the imputed data (default: ‘magic_data_allcells_sqrt.pickle’) - filtered: Set to True if data is already filtered (low quality cells and genes) (default: False)

# This command will initialize the MAGIC pipeline and generate the output directory if it does not exist
from simicpipeline import MagicPipeline
magic_pipeline = MagicPipeline(
    input_data= adata,
    project_dir='./Boiarsky_run',
    magic_output_file='magic_imputed.pickle',
    filtered=False
)

print(magic_pipeline)

Creating project directory: Boiarsky_run
MagicPipeline(
  data = AnnData object with (n_obs × n_vars) = 28550 × 22273,
  filtered = False,
  imputed = False,
  magic_data = None,
  project_dir = 'Boiarsky_run'
)

magic_pipeline.print_project_info()

Boiarsky_run/
└── magic_output/

Step 1.3: Filter Cells and Genes –> SKIP FOR BOIARSKY, already done

Remove low-quality cells and lowly-expressed genes: - min_cells_per_gene: Minimum number of cells expressing a gene (default: 10) - min_umis_per_cell: Minimum total UMI counts per cell (default: 500)

Note: If your data was already filtered you can skip this step and set the flitered argument flag to True in the previous step.

Step 1.4: Normalize Data

Perform library size normalization with scprep followed by square root transformation.

Note: this will overide adata.X with normalized data and remove adata.raw slot.

# Note this will overide adata.X with normalized data and remove adata.raw slot
magic_pipeline.normalize_data()

Normalizing data...The history saving thread hit an unexpected error (OperationalError('attempt to write a readonly database')).History will not be written to the database.

After normalization: 28550 cells x 22273 genes

MagicPipeline(
  data = AnnData object with (n_obs × n_vars) = 28550 × 22273,
  filtered = False,
  imputed = False,
  magic_data = None,
  project_dir = 'Boiarsky_run'
)

Step 1.5: Run MAGIC Imputation

Run magic imputation with defaul parameters

magic_pipeline.run_magic(
    random_state=123,
    n_jobs=-2,  # Use all but 1 CPU cores
    save_data=True
)

Warning: Data has not been filtered. Consider running filter_cells_and_genes() first.

Running MAGIC imputation...
Calculating MAGIC...
  Running MAGIC on 28550 cells and 22273 genes.
  Calculating graph and diffusion operator...
    Calculating PCA...
    Calculated PCA in 56.76 seconds.
    Calculating KNN search...
    Calculated KNN search in 9.42 seconds.
    Calculating affinities...
    Calculated affinities in 14.35 seconds.
  Calculated graph and diffusion operator in 80.59 seconds.
  Running MAGIC with `solver='exact'` on 22273-dimensional data may take a long time. Consider denoising specific genes with `genes=<list-like>` or using `solver='approximate'`.
  Calculating imputation...
  Calculated imputation in 62.68 seconds.
Calculated MAGIC in 144.20 seconds.
MAGIC imputation complete:  28550 cells x 22273 genes

Saving MAGIC-imputed data to Boiarsky_run/magic_output/magic_imputed.pickle
Saved successfully to Boiarsky_run/magic_output/magic_imputed.pickle

MagicPipeline(
  data = AnnData object with (n_obs × n_vars) = 28550 × 22273,
  filtered = False,
  imputed = True,
  magic_data = AnnData object with n_obs × n_vars = 28550 × 22273,
  project_dir = 'Boiarsky_run'
)

magic_pipeline.magic_adata.write_h5ad(magic_pipeline.magic_output_file.with_suffix('.h5ad'))

If you want to run MAGIC imputation with custom parameters you can pass them as **kwargs: - t: Number of diffusion steps (default: ‘auto’) - knn: Number of nearest neighbors (default: 5) - decay: Decay rate for kernel (default: 1) - n_jobs: Number of parallel jobs (default: -2) - genes: Genes to be returned. If None or “all genes” it returns teh entire matrix. - save_data: Whether to automatically save imputed data (default: True). If magic_output_file extension is .pickle will save it in .pickle, if h5ad, will save in adata format.

See MAGIC documentation for more parameter options.

Step 1.6: Check Pipeline Status

magic_pipeline.print_project_info(max_depth=2)

Boiarsky_run/
└── magic_output/
    ├── magic_imputed.h5ad
    └── magic_imputed.pickle

Success! MAGIC imputation is complete. The imputed data is saved in the magic_output directory.

Part 2: Experiment Setup and Gene Selection

Now we will select top variable genes and prepare input files for SimiCPipeline.

Note: if you have your data filtered, normalized and imputed with alternative methods you can start from here.

Step 2.1: Load Imputed Data

In this example we will start from the imputed AnnData object from the MAGIC pipeline. If you saved and stopped your work, you can re-load the object with the following code:

# If you saved it in h5ad format, you can load it back using:
# import simicpipeline
# imputed_data = simicpipeline.load_from_anndata('./SimiCExampleRun/magic_output/magic_imputed.h5ad')

# If you saved it in pickle format, you can load it back using:
import pickle
with open('./Boiarsky_run/magic_output/magic_imputed.pickle', 'rb') as f:
    imputed_data = pickle.load(f)
imputed_data

AnnData object with n_obs × n_vars = 28550 × 22273
    obs: 'sample_ID', 'disease_stage', 'n_genes', 'frac_mito', 'n_counts', 'cluster', 'normal_or_neoplastic'

print(f"Imputed data shape: {imputed_data.shape}")
print(imputed_data.obs.head())

Imputed data shape: (28550, 22273)
                             sample_ID disease_stage  n_genes  frac_mito  \
Cell                                                                       
AAACCTGAGGTAAACT-1-MM-5.138P      MM-5            MM     1151   0.129362   
AAACCTGGTCCCTACT-1-MM-5.138P      MM-5            MM     1260   0.108620   
AAACCTGGTTACCGAT-1-MM-5.138P      MM-5            MM     1347   0.058278   
AAACGGGAGCACACAG-1-MM-5.138P      MM-5            MM     1941   0.032493   
AAACGGGTCACAATGC-1-MM-5.138P      MM-5            MM      889   0.071260   

                              n_counts cluster normal_or_neoplastic  
Cell                                                                 
AAACCTGAGGTAAACT-1-MM-5.138P    5875.0       9           neoplastic  
AAACCTGGTCCCTACT-1-MM-5.138P    6021.0       9           neoplastic  
AAACCTGGTTACCGAT-1-MM-5.138P    8082.0       9           neoplastic  
AAACGGGAGCACACAG-1-MM-5.138P   17296.0       9           neoplastic  
AAACGGGTCACAATGC-1-MM-5.138P    3817.0       9           neoplastic  

Step 2.2: Initialize Experiment Setup

Create an experiment setup instance and directories: - input_data: Your imputed AnnData object or pandas DataFrame (cells × genes) - tf_path: Path to transcription factor (TF) list file (.csv or .txt) - project_dir: Directory where experiment files will be saved

Note: In case you do not have a TF list:

We provide a mouse TF list in the data folder that can be saved in your working data directory. - TF mouse list was downloaded in December 2024 from AnimalTFDB4 - TF human list was downloaded in February 2026 from XX

In this tutorial we are working wiht mouse data so we will use the TF list from AnimalTFDB4

# Initialize ExperimentSetup
from simicpipeline import ExperimentSetup
experiment = ExperimentSetup(
    input_data = imputed_data, 
    tf_path = "./data/JASPAR2024_Human_TFs.csv", # Should have no header, created from JASPAR2024
    project_dir='./Boiarsky_run'
)

print(f"Matrix shape: {experiment.matrix.shape}")
print(f"Number of cells: {len(experiment.cell_names)}")
print(f"Number of genes: {len(experiment.gene_names)}")
print(f"Number of TFs: {len(experiment.tf_list)}")
print(f"... Example TF names: {experiment.tf_list[0:5]}\n")
print("\n" + "="*70)
print(f"Current directory status")
print("="*70 + "\n")
experiment.print_project_info(max_depth=2)

Matrix shape: (28550, 22273)
Number of cells: 28550
Number of genes: 22273
Number of TFs: 989
... Example TF names: ['ADNP', 'ALX3', 'ALX4::TBX21', 'AR', 'ARGFX']


======================================================================
Current directory status
======================================================================

Boiarsky_run/
├── inputFiles/
├── magic_output/
│   ├── magic_imputed.h5ad
│   └── magic_imputed.pickle
└── outputSimic/
    ├── figures/
    └── matrices/

Step 2.3: Calculate MAD and Select Genes

Select top variable genes based on Median Absolute Deviation (MAD): - n_tfs: Number of top TF genes to select (default: 100) - n_targets: Number of top target genes to select (default: 1000)

Returns a tuple of (TF_list, TARGET_list)

tf_list, target_list = experiment.calculate_mad_genes(
    n_tfs=100,
    n_targets=1000
)

print(f"Selected {len(tf_list)} TFs")
print(f"Selected {len(target_list)} targets")
print(f"\nTop 10 TFs: {tf_list[:10]}")
print(f"\nTop 10 targets: {target_list[:10]}")

Removing 212 targets with MAD = 0
Selecting top 1000 targets based on MAD.
Selected 100 TFs
Selected 1000 targets

Top 10 TFs: ['JUN', 'KLF6', 'XBP1', 'JUNB', 'YBX1', 'FOS', 'KLF2', 'ATF4', 'MEF2C', 'ATF5']

Top 10 targets: ['IGKC', 'IGHA1', 'IGLC2', 'IGHG1', 'IGHG3', 'IGHG4', 'IGLC3', 'JCHAIN', 'MALAT1', 'IGHA2']

Step 2.4: Subset Data to Selected Genes

Create a subset of your data containing only the selected TFs and targets.

# Combine TF and target lists
import anndata as ad
selected_genes = tf_list + target_list

# Subset the data
if isinstance(imputed_data, ad.AnnData):
    subset_data = imputed_data[:, selected_genes].copy()
elif isinstance(imputed_data, pd.DataFrame):
    subset_data = imputed_data[selected_genes].copy()

print(f"Subset data shape: {subset_data.shape}")

Subset data shape: (28550, 1100)

Step 2.5: Save Experiment Files

Save the expression matrix and TF names in .pickle format and annotation file (optional) as .txt - run_data: ad.AnnData or pd.Dataframe with data to run in SimiC (Inputed and sliced according to experiment run) - matrix_filename: Filename to save run_data (saved with row/column headers). Can be .pickleor csv. - tf_filename: Filename for TF names list for the experiment run. Can be .pickleor csv. Even though you have a general TF_list file, this function will save the TFs selected by MAD that are found in your run_data.

• annotation:str (Optional) if run_data is ad.AnnData and annotation is in run_data.obs.columns, it will create a .csv file with the phenotype annotations needed for SimiC with cell names as index and columns category and labels.

• annotation_order(Optional) List defining the desired order of annotation categories (e.g. [‘control’, ‘treated’]). The first element maps to 0, second to 1, etc. If None, pd.factorize default order is used.

All files are saved in the inputFiles/ directory.

We recommend saving it in pickle format for fast load/dump process and save disk space.

Regarding label handling. Because it is important for SimiCPipeline run to clearly define the desired labels in the correct order, several checks and warnings will be raised if annotation_order is wrongly provided.

experiment.save_experiment_files(
    run_data = subset_data,
    matrix_filename = 'expression_matrix.pickle',
    tf_filename = 'TF_final.csv', # Will raise warning if file already exists and overwrite
    annotation = 'disease_stage', # Will rase WARNING because it is categorical, it will convert it to numerical and assign the order with the annotation_order argument
    annotation_order = ['NBM', 'SMM', 'MM'] # Will raise ERROR if missing categories or if categories in annotation_order do not match those in subset_data.obs['treatment']
)

Saved expression matrix to Boiarsky_run/inputFiles/expression_matrix.pickle
Saved 100 TFs to Boiarsky_run/inputFiles/TF_final.csv

-------

Annotation 'disease_stage' found in obs columns!
Warning: annotation is not numeric. Will convert from categorical to numeric.
Annotation order applied: {0: 'NBM', 1: 'SMM', 2: 'MM'}

Annotation distribution:
label
0     9329
1     8431
2    10790
Saved annotation to Boiarsky_run/inputFiles/disease_stage_annotation.csv

-------

Experiment files saved successfully.

-------

Be careful! Double check your labels are correct and match your cell numbers

Success! All preprocessing steps completed. Your files are ready for SimiC analysis.

Step 2.6: Verify Saved Files

Check that all files were created correctly.

experiment.print_project_info(max_depth=2)

Boiarsky_run/
├── inputFiles/
│   ├── TF_final.csv
│   ├── disease_stage_annotation.csv
│   └── expression_matrix.pickle
├── magic_output/
│   ├── magic_imputed.h5ad
│   └── magic_imputed.pickle
└── outputSimic/
    ├── figures/
    └── matrices/

Summary

This tutorial covered:

✓ Loading and filtering scRNA-seq data

✓ Running MAGIC imputation

✓ Selecting top variable genes using MAD

✓ Preparing input files for SimiC analysis with proper directory structure

Output Directory Structure

Your output directory now contains:

SimicExampleRun/
├── magic_output/
│   └── magic_imputed.pickle/.h5ad
├── inputFiles/
│   ├── expression_matrix.pickle/.csv
│   ├── TF_list.pickle
│   └── treatment_annotation.csv
└── outputSimic/
    ├── figures/
    └── matrices/

Next steps:

1. Run SimiC: Use SimiCPipeline class to run SimiC.
2. Explore results: Use SimicVisualization class to analyze GRNs and TF activities.

Check Tutorial_SimiCPipeline_full.ipynb or Tutorial_SimiCPipeline_visualization for guided info.

Final Notes

Data Format: All matrices are stored as cells × genes (rows = cells, columns = genes)

Memory Usage: MAGIC imputation can be memory-intensive for large datasets. Consider using a machine with sufficient RAM and adjusting MAGIC parameters (n_jobs, knn, t)

Please note: Although you will be able to pass custom file/direcotry paths, we highly recommend to follow the directory structure described above and follow this tutorial before running SimiC to avoid errors.