Preprocess GRN

Contents

Preprocess GRN#

Inject pySCENIC learned GRN into RegVelo

Library imports#

import numpy as np
import pandas as pd

import anndata as ad
import scanpy as sc
import scvelo as scv
from velovi import preprocess_data

from rgv_tools import DATA_DIR
from rgv_tools.preprocessing import filter_genes, set_prior_grn

/home/icb/weixu.wang/miniconda3/envs/regvelo_test/lib/python3.10/site-packages/anndata/utils.py:429: FutureWarning: Importing read_csv from `anndata` is deprecated. Import anndata.io.read_csv instead.
  warnings.warn(msg, FutureWarning)
/home/icb/weixu.wang/miniconda3/envs/regvelo_test/lib/python3.10/site-packages/anndata/utils.py:429: FutureWarning: Importing read_excel from `anndata` is deprecated. Import anndata.io.read_excel instead.
  warnings.warn(msg, FutureWarning)
/home/icb/weixu.wang/miniconda3/envs/regvelo_test/lib/python3.10/site-packages/anndata/utils.py:429: FutureWarning: Importing read_hdf from `anndata` is deprecated. Import anndata.io.read_hdf instead.
  warnings.warn(msg, FutureWarning)
/home/icb/weixu.wang/miniconda3/envs/regvelo_test/lib/python3.10/site-packages/anndata/utils.py:429: FutureWarning: Importing read_loom from `anndata` is deprecated. Import anndata.io.read_loom instead.
  warnings.warn(msg, FutureWarning)
/home/icb/weixu.wang/miniconda3/envs/regvelo_test/lib/python3.10/site-packages/anndata/utils.py:429: FutureWarning: Importing read_mtx from `anndata` is deprecated. Import anndata.io.read_mtx instead.
  warnings.warn(msg, FutureWarning)
/home/icb/weixu.wang/miniconda3/envs/regvelo_test/lib/python3.10/site-packages/anndata/utils.py:429: FutureWarning: Importing read_text from `anndata` is deprecated. Import anndata.io.read_text instead.
  warnings.warn(msg, FutureWarning)
/home/icb/weixu.wang/miniconda3/envs/regvelo_test/lib/python3.10/site-packages/anndata/utils.py:429: FutureWarning: Importing read_umi_tools from `anndata` is deprecated. Import anndata.io.read_umi_tools instead.
  warnings.warn(msg, FutureWarning)

Constants#

DATASET = "mouse_neural_crest"

SAVE_DATA = True
if SAVE_DATA:
    (DATA_DIR / DATASET / "processed").mkdir(parents=True, exist_ok=True)

Data loading#

TF_all = pd.read_csv("/home/icb/weixu.wang/regulatory_velo/pancreas_dataset/allTFs_mm.txt", header=None)
ldata = ad.io.read_h5ad(DATA_DIR / DATASET / "raw" / "GSE201257_adata_velo_raw.h5ad")

Preprocess for each GRN#

for i in range(1, 5):
    print("loading")
    reg = pd.read_csv(DATA_DIR / DATASET / "processed" / f"regulon_mat_stage_{i}_all_regulons.csv", index_col=0)

    reg.index = reg.index.str.extract(r"(\w+)")[0]
    reg = reg.groupby(reg.index).sum()
    reg[reg != 0] = 1

    TF = np.unique([x.split("(")[0] for x in reg.index.tolist()])
    genes = np.unique(TF.tolist() + reg.columns.tolist())

    GRN = pd.DataFrame(0, index=genes, columns=genes)
    GRN.loc[TF, reg.columns.tolist()] = np.array(reg)

    mask = (GRN.sum(0) != 0) | (GRN.sum(1) != 0)
    GRN = GRN.loc[mask, mask].copy()

    if SAVE_DATA:
        GRN.to_parquet(DATA_DIR / DATASET / "processed" / f"regulon_mat_stage_{i}_processed_all_regulons.parquet")
    print("Done! processed GRN with " + str(reg.shape[0]) + " TF and " + str(reg.shape[1]) + " targets")

loading
Done! processed GRN with 738 TF and 19054 targets
loading
Done! processed GRN with 722 TF and 20164 targets
loading
Done! processed GRN with 724 TF and 20843 targets
loading
Done! processed GRN with 730 TF and 20843 targets

Preprocess anndata#

for i in range(1, 5):
    print("loading...")
    adata = scv.read(DATA_DIR / DATASET / "processed" / f"adata_stage{i}_processed.h5ad")
    adata = scv.utils.merge(adata, ldata)
    adata.obsp = None
    del adata.uns["neighbors"]

    print(adata)
    adata.var["TF"] = [i in TF_all.iloc[:, 0].tolist() for i in adata.var_names.tolist()]

    ## select highly variable genes
    scv.pp.filter_genes(adata, min_shared_counts=20)
    scv.pp.normalize_per_cell(adata)

    ## In top 2000 TF is too few, we select all TFs in top 3000
    sc.pp.highly_variable_genes(adata, n_top_genes=2000)
    adata = adata[:, adata.var["highly_variable"]].copy()

    sc.pp.neighbors(adata, n_neighbors=30, n_pcs=30)
    scv.pp.moments(adata, n_pcs=None, n_neighbors=None)

    ## load GRN
    GRN = pd.read_parquet(DATA_DIR / DATASET / "processed" / f"regulon_mat_stage_{i}_processed_all_regulons.parquet")

    adata = set_prior_grn(adata, GRN.T)
    velocity_genes = preprocess_data(adata.copy()).var_names.tolist()
    TF = adata.var_names[adata.uns["skeleton"].sum(1) != 0]
    var_mask = np.union1d(TF, velocity_genes)

    adata = adata[:, var_mask].copy()
    adata = filter_genes(adata)
    adata = preprocess_data(adata, filter_on_r2=False)

    adata.var["velocity_genes"] = adata.var_names.isin(velocity_genes)
    adata.var["TF"] = adata.var_names.isin(TF)

    print(adata)

    if SAVE_DATA:
        adata.write_h5ad(DATA_DIR / DATASET / "processed" / f"adata_stage{i}_processed_velo_all_regulons.h5ad")

loading...
AnnData object with n_obs × n_vars = 3926 × 20559
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced'
Filtered out 5602 genes that are detected 20 counts (shared).
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set `enforce=True`.
Normalized count data: spliced, unspliced.

2025-08-19 15:54:36.626136: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:477] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
E0000 00:00:1755611677.677800 2423155 cuda_dnn.cc:8310] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
E0000 00:00:1755611678.146442 2423155 cuda_blas.cc:1418] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered

computing moments based on connectivities
    finished (0:00:00) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
computing velocities
    finished (0:00:00) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)
Number of genes: 1205
Number of genes: 1193
AnnData object with n_obs × n_vars = 3926 × 1193
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'TF', 'means', 'dispersions', 'dispersions_norm', 'velocity_genes'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap', 'hvg', 'neighbors', 'regulators', 'targets', 'skeleton', 'network'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced', 'Ms', 'Mu'
    obsp: 'distances', 'connectivities'
loading...
AnnData object with n_obs × n_vars = 5139 × 20559
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced'
Filtered out 4339 genes that are detected 20 counts (shared).
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set `enforce=True`.
Normalized count data: spliced, unspliced.
computing moments based on connectivities
    finished (0:00:00) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
computing velocities
    finished (0:00:00) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)
Number of genes: 1160
Number of genes: 1124
AnnData object with n_obs × n_vars = 5139 × 1124
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'TF', 'means', 'dispersions', 'dispersions_norm', 'velocity_genes'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap', 'hvg', 'neighbors', 'regulators', 'targets', 'skeleton', 'network'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced', 'Ms', 'Mu'
    obsp: 'distances', 'connectivities'
loading...
AnnData object with n_obs × n_vars = 6273 × 20559
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced'
Filtered out 4178 genes that are detected 20 counts (shared).
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set `enforce=True`.
Normalized count data: spliced, unspliced.
computing moments based on connectivities
    finished (0:00:00) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
computing velocities
    finished (0:00:00) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)
Number of genes: 1159
Number of genes: 1079
Number of genes: 1065
AnnData object with n_obs × n_vars = 6273 × 1065
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'TF', 'means', 'dispersions', 'dispersions_norm', 'velocity_genes'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap', 'hvg', 'neighbors', 'regulators', 'targets', 'skeleton', 'network'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced', 'Ms', 'Mu'
    obsp: 'distances', 'connectivities'
loading...
AnnData object with n_obs × n_vars = 8821 × 20559
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced'
Filtered out 3760 genes that are detected 20 counts (shared).
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set `enforce=True`.
Normalized count data: spliced, unspliced.
computing moments based on connectivities
    finished (0:00:00) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
computing velocities
    finished (0:00:00) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)
Number of genes: 1176
Number of genes: 1150
AnnData object with n_obs × n_vars = 8821 × 1150
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'all_states', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'TF', 'means', 'dispersions', 'dispersions_norm', 'velocity_genes'
    uns: 'all_states_colors', 'assignments_colors', 'devtime_colors', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'paga', 'umap', 'hvg', 'neighbors', 'regulators', 'targets', 'skeleton', 'network'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    layers: 'GEX', 'palantir_imp', 'scaled', 'ambiguous', 'matrix', 'spanning', 'spliced', 'unspliced', 'Ms', 'Mu'
    obsp: 'distances', 'connectivities'