Dynamo-based perturbation analysis

Notebooks for predicting TF perturbation effects using dynamo.

Library imports

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import mplscience
import seaborn as sns

import dynamo as dyn
import scanpy as sc
from dynamo.preprocessing import Preprocessor

from rgv_tools import DATA_DIR, FIG_DIR
from rgv_tools.perturbation import (
    abundance_test,
    get_list_name,
    Multiple_TFScanning_KO_dyn,
    Multiple_TFScanning_perturbation_dyn,
    split_elements,
    TFScanning_KO_dyn,
    TFScanning_perturbation_dyn,
)

/home/icb/weixu.wang/miniconda3/envs/dynamo/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/dynamo/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/dynamo/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/dynamo/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/dynamo/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/dynamo/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/dynamo/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 = "zebrafish"

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

SAVE_FIGURES = False
if SAVE_FIGURES:
    (FIG_DIR / DATASET).mkdir(parents=True, exist_ok=True)

genes = ["nr2f5", "sox9b", "twist1b", "ets1"]

TERMINAL_STATES = [
    "mNC_head_mesenchymal",
    "mNC_arch2",
    "mNC_hox34",
    "Pigment",
]

single_ko = ["elk3", "erf", "fli1a", "mitfa", "nr2f5", "rarga", "rxraa", "smarcc1a", "tfec", "nr2f2"]
multiple_ko = ["fli1a_elk3", "mitfa_tfec", "tfec_mitfa_bhlhe40", "fli1a_erf_erfl3", "erf_erfl3"]

Data loading

adata = sc.read_h5ad(DATA_DIR / DATASET / "processed" / "adata_preprocessed.h5ad")

adata.X = adata.layers["matrix"].copy()

Processing by dynamo

preprocessor = Preprocessor()
preprocessor.preprocess_adata(adata, recipe="monocle")

|-----> Running monocle preprocessing pipeline...
|-----------> filtered out 0 outlier cells
|-----------> filtered out 21 outlier genes
|-----> PCA dimension reduction
|-----> <insert> X_pca to obsm in AnnData Object.
|-----> [Preprocessor-monocle] completed [1.1507s]

dyn.tl.dynamics(adata)

|-----> dynamics_del_2nd_moments_key is None. Using default value from DynamoAdataConfig: dynamics_del_2nd_moments_key=False
|-----------> removing existing M layers:[]...
|-----------> making adata smooth...
|-----> calculating first/second moments...
|-----? layer X_velocity is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----? layer X_variance_velocity is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----? layer X_Mu is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----? layer X_Ms is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----? layer X_Mu is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----? layer X_Ms is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----? layer X_Mu is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----? layer X_Ms is not in any of the (['X_spliced', 'X_unspliced'], ['X_new', 'X_total'], ['X_uu', 'X_ul', 'X_su', 'X_sl']) groups, skipping...
|-----> [moments calculation] completed [11.9131s]

estimating gamma: 100%|██████████| 988/988 [00:10<00:00, 96.84it/s] 

AnnData object with n_obs × n_vars = 697 × 988
    obs: 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts', 'cell_type', 'stage', 'cell_type2', 'macrostates', 'nGenes', 'nCounts', 'pMito', 'pass_basic_filter', 'unspliced_Size_Factor', 'initial_unspliced_cell_size', 'Size_Factor', 'initial_cell_size', 'spliced_Size_Factor', 'initial_spliced_cell_size', 'variance_velocity_Size_Factor', 'initial_variance_velocity_cell_size', 'Ms_Size_Factor', 'initial_Ms_cell_size', 'velocity_Size_Factor', 'initial_velocity_cell_size', 'Mu_Size_Factor', 'initial_Mu_cell_size', 'ntr'
    var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'gene_count_corr', 'is_tf', 'velocity_gamma', 'velocity_qreg_ratio', 'velocity_r2', 'velocity_genes', 'TF', 'nCells', 'nCounts', 'pass_basic_filter', 'frac', 'use_for_pca', 'ntr', 'beta', 'gamma', 'half_life', 'alpha_b', 'alpha_r2', 'gamma_b', 'gamma_r2', 'gamma_logLL', 'delta_b', 'delta_r2', 'bs', 'bf', 'uu0', 'ul0', 'su0', 'sl0', 'U0', 'S0', 'total0', 'use_for_dynamics'
    uns: 'cell_type_colors', 'macrostates_colors', 'neighbors', 'network', 'regulators', 'skeleton', 'targets', 'velocity_params', 'pp', 'feature_selection', 'PCs', 'explained_variance_ratio_', 'pca_mean', 'dynamics'
    obsm: 'X_pca', 'X_umap'
    layers: 'Ms', 'Mu', 'ambiguous', 'matrix', 'spliced', 'unspliced', 'variance_velocity', 'velocity', 'X_spliced', 'X_variance_velocity', 'X_unspliced', 'X_velocity', 'X_Mu', 'X_Ms', 'M_u', 'M_uu', 'M_s', 'M_us', 'M_ss', 'velocity_S'
    obsp: 'connectivities', 'distances', 'moments_con'

dyn.tl.reduceDimension(adata)

|-----> retrieve data for non-linear dimension reduction...
|-----? adata already have basis umap. dimension reduction umap will be skipped! 
set enforce=True to re-performing dimension reduction.
|-----> [UMAP] completed [0.0022s]

dyn.tl.cell_velocities(adata, basis="pca")

|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.5832s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.3028s]

AnnData object with n_obs × n_vars = 697 × 988
    obs: 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts', 'cell_type', 'stage', 'cell_type2', 'macrostates', 'nGenes', 'nCounts', 'pMito', 'pass_basic_filter', 'unspliced_Size_Factor', 'initial_unspliced_cell_size', 'Size_Factor', 'initial_cell_size', 'spliced_Size_Factor', 'initial_spliced_cell_size', 'variance_velocity_Size_Factor', 'initial_variance_velocity_cell_size', 'Ms_Size_Factor', 'initial_Ms_cell_size', 'velocity_Size_Factor', 'initial_velocity_cell_size', 'Mu_Size_Factor', 'initial_Mu_cell_size', 'ntr'
    var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'gene_count_corr', 'is_tf', 'velocity_gamma', 'velocity_qreg_ratio', 'velocity_r2', 'velocity_genes', 'TF', 'nCells', 'nCounts', 'pass_basic_filter', 'frac', 'use_for_pca', 'ntr', 'beta', 'gamma', 'half_life', 'alpha_b', 'alpha_r2', 'gamma_b', 'gamma_r2', 'gamma_logLL', 'delta_b', 'delta_r2', 'bs', 'bf', 'uu0', 'ul0', 'su0', 'sl0', 'U0', 'S0', 'total0', 'use_for_dynamics', 'use_for_transition'
    uns: 'cell_type_colors', 'macrostates_colors', 'neighbors', 'network', 'regulators', 'skeleton', 'targets', 'velocity_params', 'pp', 'feature_selection', 'PCs', 'explained_variance_ratio_', 'pca_mean', 'dynamics', 'grid_velocity_pca'
    obsm: 'X_pca', 'X_umap', 'velocity_pca'
    layers: 'Ms', 'Mu', 'ambiguous', 'matrix', 'spliced', 'unspliced', 'variance_velocity', 'velocity', 'X_spliced', 'X_variance_velocity', 'X_unspliced', 'X_velocity', 'X_Mu', 'X_Ms', 'M_u', 'M_uu', 'M_s', 'M_us', 'M_ss', 'velocity_S'
    obsp: 'connectivities', 'distances', 'moments_con', 'pearson_transition_matrix'

dyn.vf.VectorField(adata, basis="pca")

|-----> VectorField reconstruction begins...
|-----> Retrieve X and V based on basis: PCA. 
        Vector field will be learned in the PCA space.
|-----> Learning vector field with method: sparsevfc.
|-----> [SparseVFC] begins...
|-----> Sampling control points based on data velocity magnitude...
|-----> [SparseVFC] completed [0.1569s]
|-----> [VectorField] completed [0.2136s]

adata_perturb = adata.copy()

del adata.uns["cell_type_colors"]

Compare predict fate probability changes

fate_prob = {}
fate_prob_perturb = {}

for g in genes:
    fb, fb_perturb = TFScanning_KO_dyn(adata, 8, "cell_type", TERMINAL_STATES, [g], fate_prob_return=True)
    fate_prob[g] = fb
    fate_prob_perturb[g] = fb_perturb

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:13,722 - INFO - Using pre-computed Schur decomposition

WARNING: Unable to import petsc4py. For installation, please refer to: https://petsc4py.readthedocs.io/en/stable/install.html.
Defaulting to `'gmres'` solver.

|-----> In silico knockout nr2f5
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.3101s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2672s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:18,010 - INFO - Using pre-computed Schur decomposition

Done nr2f5

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:20,505 - INFO - Using pre-computed Schur decomposition

|-----> In silico knockout sox9b
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2859s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2961s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:23,720 - INFO - Using pre-computed Schur decomposition

Done sox9b

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:27,163 - INFO - Using pre-computed Schur decomposition

|-----> In silico knockout twist1b
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2847s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2827s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:31,388 - INFO - Using pre-computed Schur decomposition

Done twist1b

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:34,859 - INFO - Using pre-computed Schur decomposition

|-----> In silico knockout ets1
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2764s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2755s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:38,796 - INFO - Using pre-computed Schur decomposition

Done ets1

Visualize dynamo perturbation effects

for g in genes:
    data = abundance_test(fate_prob[g], fate_prob_perturb[g])
    data = pd.DataFrame(
        {
            "Score": data.iloc[:, 0].tolist(),
            "p-value": data.iloc[:, 1].tolist(),
            "Terminal state": data.index.tolist(),
            "TF": [g] * (data.shape[0]),
        }
    )

    final_df = data.copy()
    final_df["Score"] = 0.5 - final_df["Score"]

    color_label = "cell_type"
    df = pd.DataFrame(final_df["Score"])
    df.columns = ["coefficient"]
    df["Cell type"] = final_df["Terminal state"]
    order = df["Cell type"].tolist()

    palette = dict(zip(adata.obs[color_label].cat.categories, adata_perturb.uns[f"{color_label}_colors"]))
    subset_palette = {name: color for name, color in palette.items() if name in TERMINAL_STATES}

    with mplscience.style_context():
        sns.set_style(style="whitegrid")
        fig, ax = plt.subplots(figsize=(2, 2))
        sns.barplot(
            data=df,
            y="coefficient",
            x="Cell type",
            palette=subset_palette,
            order=order,
            ax=ax,
        )
        ax.tick_params(axis="x", rotation=90)
        plt.title("$\\mathit{" + g + "}$ regulon knock out simulation")

        if SAVE_FIGURES:
            plt.savefig(
                FIG_DIR / DATASET / f"{g}_perturbation_simulation_dynamo.svg",
                format="svg",
                transparent=True,
                bbox_inches="tight",
            )
        # Show the plot
        plt.show()

Perturbation prediction

Single gene knockout

single_ko = set(single_ko).intersection(adata.var_names)
single_ko = list(single_ko)

## Dynamo (KO)
ko_dyn = TFScanning_KO_dyn(adata, 8, "cell_type", TERMINAL_STATES, single_ko)

## Dynamo (perturbation)
perturbation_dyn = TFScanning_perturbation_dyn(adata, 8, "cell_type", TERMINAL_STATES, single_ko)

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:44,048 - INFO - Using pre-computed Schur decomposition

|-----> In silico knockout nr2f5
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2852s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2654s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:48,132 - INFO - Using pre-computed Schur decomposition

Done nr2f5
|-----> In silico knockout nr2f2
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.1915s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.1805s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:51,004 - INFO - Using pre-computed Schur decomposition

Done nr2f2
|-----> In silico knockout elk3
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2922s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2787s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:55,360 - INFO - Using pre-computed Schur decomposition

Done elk3
|-----> In silico knockout mitfa
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.3077s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2411s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:02:58,967 - INFO - Using pre-computed Schur decomposition

Done mitfa
|-----> In silico knockout fli1a
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2957s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2917s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:02,114 - INFO - Using pre-computed Schur decomposition

Done fli1a
|-----> In silico knockout rxraa
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2815s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2870s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:06,770 - INFO - Using pre-computed Schur decomposition

Done rxraa
|-----> In silico knockout erf
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2904s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2832s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:11,135 - INFO - Using pre-computed Schur decomposition

Done erf
|-----> In silico knockout rarga
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2852s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2603s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:15,267 - INFO - Using pre-computed Schur decomposition

Done rarga
|-----> In silico knockout tfec
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2959s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2998s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:18,906 - INFO - Using pre-computed Schur decomposition

Done tfec
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....

calculating Jacobian for each cell: 100%|██████████| 697/697 [00:00<00:00, 109742.48it/s]

|-----> project the pca perturbation vector to low dimensional space....

|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.3244s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done nr2f5
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2364s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done nr2f2
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.3204s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done elk3
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.3849s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done mitfa
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2955s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done fli1a
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2752s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done rxraa
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2678s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done erf
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.3549s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done rarga
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2676s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='pca_perturbation') to visualize the perturbation vector
Done tfec

## Perform KO screening using function based perturbation
coef_KO = pd.DataFrame(np.array(ko_dyn["coefficient"]))
coef_KO.index = ko_dyn["TF"]
coef_KO.columns = get_list_name(ko_dyn["coefficient"][0])

## Perform perturbation screening using gene expression based perturbation
coef_perturb = pd.DataFrame(np.array(perturbation_dyn["coefficient"]))
coef_perturb.index = perturbation_dyn["TF"]
coef_perturb.columns = get_list_name(perturbation_dyn["coefficient"][0])

Multiple gene knock-out prediction

multiple_ko_list = split_elements(multiple_ko)

## Dynamo (KO)
ko_dyn = Multiple_TFScanning_KO_dyn(adata, 9, "cell_type", TERMINAL_STATES, multiple_ko_list)

## Dynamo (perturbation)
perturbation_ko = Multiple_TFScanning_perturbation_dyn(adata, 9, "cell_type", TERMINAL_STATES, multiple_ko_list)

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:33,803 - INFO - Using pre-computed Schur decomposition

|-----> In silico knockout ['fli1a', 'elk3']
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.3093s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2854s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:38,437 - INFO - Using pre-computed Schur decomposition

['fli1a', 'elk3']
Done fli1a_elk3
|-----> In silico knockout ['mitfa', 'tfec']
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.4391s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.4197s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:44,295 - INFO - Using pre-computed Schur decomposition

['mitfa', 'tfec']
Done mitfa_tfec
|-----> In silico knockout ['tfec', 'mitfa', 'bhlhe40']
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2762s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2686s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:48,604 - INFO - Using pre-computed Schur decomposition

['tfec', 'mitfa', 'bhlhe40']
Done tfec_mitfa_bhlhe40
|-----> In silico knockout ['fli1a', 'erf', 'erfl3']
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.3052s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2681s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:52,901 - INFO - Using pre-computed Schur decomposition

['fli1a', 'erf', 'erfl3']
Done fli1a_erf_erfl3
|-----> In silico knockout ['erf', 'erfl3']
|-----> Project the high dimensional vector field after KO to umap.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%|-----> [calculating transition matrix via pearson kernel with sqrt transform.] completed [0.2932s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2666s]

WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying

2024-11-24 16:03:57,243 - INFO - Using pre-computed Schur decomposition

['erf', 'erfl3']
Done erf_erfl3
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2709s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='umap_perturbation') to visualize the perturbation vector
Done fli1a_elk3
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2763s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='umap_perturbation') to visualize the perturbation vector
Done mitfa_tfec
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2808s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='umap_perturbation') to visualize the perturbation vector
Done tfec_mitfa_bhlhe40
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2947s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='umap_perturbation') to visualize the perturbation vector
Done fli1a_erf_erfl3
|-----> In silico perturbation of single-cells and prediction of cell fate after perturbation...
|-----> Retrive X_pca, PCs, pca_mean...
|-----> Calculate perturbation effect matrix via \delta Y = J \dot \delta X....
|-----> project the pca perturbation vector to low dimensional space....
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%|-----> [projecting velocity vector to low dimensional embedding] completed [0.2571s]
|-----> you can use dyn.pl.streamline_plot(adata, basis='umap_perturbation') to visualize the perturbation vector
Done erf_erfl3

## Perform KO screening using function based perturbation
coef_KO_multiple = pd.DataFrame(np.array(ko_dyn["coefficient"]))
coef_KO_multiple.index = ko_dyn["TF"]
coef_KO_multiple.columns = get_list_name(ko_dyn["coefficient"][0])

## Perform perturbation screening using gene expression based perturbation
coef_perturb_multiple = pd.DataFrame(np.array(perturbation_dyn["coefficient"]))
coef_perturb_multiple.index = perturbation_dyn["TF"]
coef_perturb_multiple.columns = get_list_name(perturbation_dyn["coefficient"][0])

Save dataset

if SAVE_DATA:
    coef_KO.to_csv(DATA_DIR / DATASET / "results" / "dynamo_KO_single.csv")
    coef_KO_multiple.to_csv(DATA_DIR / DATASET / "results" / "dynamo_KO_multiple.csv")

    coef_perturb.to_csv(DATA_DIR / DATASET / "results" / "dynamo_perturb_single.csv")
    coef_perturb_multiple.to_csv(DATA_DIR / DATASET / "results" / "dynamo_perturb_multiple.csv")