{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
" # Lineage Driver Discovery via Pairwise KO Analysis\n",
"\n",
"\n",
"\n",
" This notebook discovers gene pairs that drive erythroid vs myeloid differentiation when knocked out,\n",
"\n",
" using the Paul et al. 2015 hematopoiesis dataset processed through scHopfield.\n",
"\n",
"\n",
"\n",
" **Pipeline overview:**\n",
"\n",
" 1. Candidate TF prioritization from W-matrix, centrality, and energy correlations\n",
"\n",
" 2. Single-gene KO baseline to rank candidates\n",
"\n",
" 3. Pairwise KO screen across top candidates\n",
"\n",
" 4. Lineage-biasing score and synergy analysis\n",
"\n",
" 5. Visualization of top pairs\n",
"\n",
"\n",
"\n",
" **Biological validation anchors:**\n",
"\n",
" - Gata1 KO → erythroid block\n",
"\n",
" - Spi1 (PU.1) KO → myeloid/GMP block\n",
"\n",
" - Klf1 KO → erythroid differentiation block\n",
"\n",
" - Gata1 + Spi1 → known antagonist pair"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" ## 0. Setup & Data Loading"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Package vs Notebook: Methods Classification\n",
"\n",
"Which analytical methods in this notebook are **generic enough to move into the `scHopfield` package** vs **study-specific and should remain here**.\n",
"\n",
"---\n",
"\n",
"### MOVE TO PACKAGE (`sch`)\n",
"\n",
"| Cell / Section | Method | Suggested API |\n",
"|----------------|--------|---------------|\n",
"| §2 — Single-gene KO screen (cell 24) | ODE-based per-gene lineage bias screen | `sch.dyn.run_ko_screen` — **already in package** |\n",
"| §3 — Pairwise KO screen (cell 34) | ODE-based double-KO lineage bias for arbitrary gene pairs | `sch.dyn.run_pairwise_ko_screen` — **already in package** |\n",
"| §1 — TF scoring (cell 18) | Composite rank-sum score (energy-gene correlation, out-degree centrality, velocity alignment) | `sch.tl.score_driver_tfs` — **already in package** |\n",
"| §4 — Epistasis metrics (cell 39) | Cancellation error = double_bias - (bias_A + bias_B); dominant epistasis = max-magnitude per-lineage deviation from best single agent | `sch.dyn.compute_epistasis(pair_df)` — **should be added to package**; currently computed inline |\n",
"| §2 — WT Hopfield velocity (cell 23) | Pre-computation of WT flow in embedding, reused for all KO comparisons | `sch.tl.compute_wt_flow(adata)` — already implicitly in package; expose as explicit checkpoint method |\n",
"\n",
"---\n",
"\n",
"### KEEP IN NOTEBOOK (study-specific)\n",
"\n",
"| Cell / Section | Method | Why notebook-only |\n",
"|----------------|--------|-------------------|\n",
"| §1 — Pareto multi-round TF selection (cell 21) | Iterative multi-objective Pareto ranking with `NUM_PARETO_ROUNDS=9`, erythroid-exclusive/myeloid-exclusive logic, cross-round candidate persistence | Specific algorithm and round count for this study; not generalisable as-is |\n",
"| §2 — Top-N selection for pairwise screen (cell 26) | Choosing top-5 erythroid and top-5 myeloid drivers from single-KO results | Cutoff choice (N=5) and cross/within lineage pair construction are study-specific |\n",
"| §4 — Scatter plots and staircase heatmap (cells 40-41) | Joint visualisation of epistasis for all 45 pairs | Specific to 45-pair screen structure and Paul 2015 cluster layout |\n",
"| §5a — UMAP flow plots for top pairs (cells 44-45) | Per-pair KO flow re-simulation and embedding visualisation | Specific gene pairs and embedding; not generalisable |\n",
"| §5b — Energy change per cluster (cells 47-48) | ΔE per cluster for top double-KO pairs | Specific pair set; general ΔE computation is in `sch.tl.compute_energies` |\n",
"| §5c — Biological validation table (cell 49) | Mapping known biology (literature) pairs to model predictions | Entirely study-specific |\n",
"| All figure generation / assembly | Matplotlib figure code, panel layout, caption writing | Study-specific |\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Device: mps\n",
"Erythroid clusters : ['1Ery', '2Ery', '3Ery', '4Ery', '5Ery', '6Ery', '7MEP']\n",
"Myeloid clusters : ['9GMP', '10GMP', '11DC', '12Baso', '13Baso', '14Mo', '15Mo', '16Neu', '17Neu', '18Eos']\n",
"Progenitor clusters : ['7MEP', '9GMP', '10GMP']\n"
]
}
],
"source": [
"import warnings\n",
"import itertools\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.colors as mcolors\n",
"import numpy as np\n",
"import pandas as pd\n",
"import scanpy as sc\n",
"import scvelo as scv\n",
"import seaborn as sns\n",
"import torch\n",
"\n",
"import celloracle as co\n",
"import scHopfield as sch\n",
"\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"# ── Load or run the pipeline ─────────────────────────────────────────────────\n",
"LOAD_DATA = True\n",
"# ── Analysis parameters ──────────────────────────────────────────────────────\n",
"CLUSTER_KEY = 'paul15_clusters'\n",
"SPLICED_KEY = 'Ms'\n",
"DEGRADATION_KEY = 'gamma'\n",
"BASIS = 'draw_graph_fa'\n",
"VELOCITY_KEY = 'velocity_S'\n",
"VELOCITY_SCALE = 500.0\n",
"SCAFFOLD_REG = 1e-1\n",
"N_EPOCHS = 1000\n",
"BATCH_SIZE = 128\n",
"DEVICE = 'cuda' if torch.cuda.is_available() else 'mps' if torch.backends.mps.is_available() else 'cpu'\n",
"\n",
"\n",
"# ── Lineage group definitions ────────────────────────────────────────────────\n",
"ERYTHROID = ['1Ery', '2Ery', '3Ery', '4Ery', '5Ery', '6Ery', '7MEP']\n",
"MYELOID = ['9GMP', '10GMP', '11DC', '12Baso', '13Baso', '14Mo', '15Mo', '16Neu', '17Neu', '18Eos']\n",
"PROGENITOR = ['7MEP', '9GMP', '10GMP'] # bifurcation-point cells\n",
"\n",
"# ── Paul15 cluster display order ─────────────────────────────────────────────\n",
"CLUSTER_ORDER = [\n",
" '1Ery', '2Ery', '3Ery', '4Ery', '5Ery', '6Ery', '7MEP',\n",
" '8Mk', '9GMP', '10GMP', '11DC', '12Baso', '13Baso',\n",
" '14Mo', '15Mo', '16Neu', '17Neu', '18Eos', '19Lymph',\n",
"]\n",
"\n",
"print(f\"Device: {DEVICE}\")\n",
"print(f\"Erythroid clusters : {ERYTHROID}\")\n",
"print(f\"Myeloid clusters : {MYELOID}\")\n",
"print(f\"Progenitor clusters : {PROGENITOR}\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" ### Checkpoint paths\n",
"\n",
"\n",
"\n",
" Three checkpoints are saved during the run.\n",
"\n",
" To **resume** from any checkpoint, skip the preceding computation cells and run the corresponding *Resume* cell instead.\n",
"\n",
"\n",
"\n",
" | Checkpoint | Saves | Covers |\n",
"\n",
" |---|---|---|\n",
"\n",
" | `adata_schopfield.h5ad` | base AnnData after GRN inference | ~30–60 min |\n",
"\n",
" | `single_ko_*.csv` | single-KO summary tables | ~15 min |\n",
"\n",
" | `pair_ko_*.csv` | pairwise-KO summary tables | ~45 min |"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib\n",
"matplotlib.rcParams._get = matplotlib.rcParams.get\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Checkpoint directory: /Users/bernaljp/Documents/SCHData/checkpoints/06_lineage_drivers\n"
]
}
],
"source": [
"import os, json\n",
"\n",
"SAVE_DIR = '/Users/bernaljp/Documents/SCHData/checkpoints/06_lineage_drivers'\n",
"os.makedirs(SAVE_DIR, exist_ok=True)\n",
"\n",
"# ── File paths ────────────────────────────────────────────────────────────────\n",
"ADATA_PATH = f'{SAVE_DIR}/adata_schopfield.h5ad'\n",
"MODEL_PATH = f'{SAVE_DIR}/model.h5sch'\n",
"TF_NAMES_PATH = f'{SAVE_DIR}/tf_names.json'\n",
"COLORS_PATH = f'{SAVE_DIR}/colors.csv'\n",
"\n",
"SINGLE_KO_BIAS_PATH = f'{SAVE_DIR}/single_ko_bias.csv'\n",
"SINGLE_KO_PER_CL_PATH = f'{SAVE_DIR}/single_ko_per_cluster.csv'\n",
"CANDIDATES_PATH = f'{SAVE_DIR}/candidates.json'\n",
"\n",
"PAIR_KO_BIAS_PATH = f'{SAVE_DIR}/pair_ko_bias.csv'\n",
"PAIR_KO_PER_CL_PATH = f'{SAVE_DIR}/pair_ko_per_cluster.csv'\n",
"\n",
"print(f\"Checkpoint directory: {os.path.abspath(SAVE_DIR)}\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" ### Rebuild the scHopfield pipeline (same as notebook 05)\n",
"\n",
"\n",
"\n",
" We re-run the full setup: load the CellOracle tutorial data, estimate velocities,\n",
"\n",
" fit sigmoids, build the GRN scaffold, and infer cluster-specific interaction matrices.\n",
"\n",
" After this cell block the `adata` object carries `varp['W_{cluster}']` and all\n",
"\n",
" downstream fields required for candidate scoring."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Loaded: 2671 cells × 1999 genes\n",
"Clusters: ['10GMP', '11DC', '12Baso', '13Baso', '14Mo', '15Mo', '16Neu', '17Neu', '18Eos', '19Lymph', '1Ery', '2Ery', '3Ery', '4Ery', '5Ery', '6Ery', '7MEP', '8Mk', '9GMP']\n"
]
}
],
"source": [
"# ── Load data ─────────────────────────────────────────────────────────────────\n",
"oracle_demo = co.data.load_tutorial_oracle_object()\n",
"adata = oracle_demo.adata.copy()\n",
"adata.var['scHopfield_used'] = True\n",
"\n",
"print(f\"Loaded: {adata.n_obs} cells × {adata.n_vars} genes\")\n",
"print(f\"Clusters: {sorted(adata.obs[CLUSTER_KEY].unique())}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"computing neighbors\n",
" finished (0:00:02) --> added \n",
" 'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)\n",
"computing moments based on connectivities\n",
" finished (0:00:00) --> added \n",
" 'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)\n",
"computing velocity graph (using 14/14 cores)\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "08c5b229175a4fb79f17204fa6976ba1",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/2671 [00:00 added \n",
" 'velocity_S_graph', sparse matrix with cosine correlations (adata.uns)\n",
"computing velocity embedding\n",
" finished (0:00:00) --> added\n",
" 'velocity_S_draw_graph_fa', embedded velocity vectors (adata.obsm)\n",
"Gamma range: [0.100, 0.100]\n"
]
}
],
"source": [
"# ── Velocity estimation ───────────────────────────────────────────────────────\n",
"adata.layers['spliced'] = adata.layers['normalized_count']\n",
"adata.layers['unspliced'] = adata.layers['normalized_count']\n",
"scv.pp.moments(adata, n_pcs=30, n_neighbors=30)\n",
"_ = adata.layers.pop('unspliced')\n",
"\n",
"sch.pp.estimate_velocity_from_pseudotime(\n",
" adata,\n",
" pseudotime_key='Pseudotime',\n",
" spliced_key=SPLICED_KEY,\n",
" connectivity_key='connectivities',\n",
" scale=VELOCITY_SCALE,\n",
" store_key=VELOCITY_KEY,\n",
")\n",
"\n",
"scv.tl.velocity_graph(adata, vkey=VELOCITY_KEY, xkey=SPLICED_KEY, n_jobs=-1)\n",
"scv.tl.velocity_embedding(adata, basis=BASIS, vkey=VELOCITY_KEY)\n",
"adata.obsm[f'velocity_{BASIS}'] = adata.obsm[f'{VELOCITY_KEY}_{BASIS}']\n",
"\n",
"# Degradation rates\n",
"expression = adata.layers[SPLICED_KEY].copy()\n",
"velocities = adata.layers[VELOCITY_KEY]\n",
"mean_expr = np.abs(expression).mean(axis=0) + 1e-6\n",
"mean_vel = np.abs(velocities).mean(axis=0)\n",
"gamma = np.clip(mean_vel / mean_expr, 0.1, 10.0)\n",
"adata.var[DEGRADATION_KEY] = gamma\n",
"\n",
"print(f\"Gamma range: [{gamma.min():.3f}, {gamma.max():.3f}]\")\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"if not LOAD_DATA:\n",
" # ── Sigmoid fitting ────────────────────────────────────────────────────────────\n",
" adata.var['scHopfield_used'] = True\n",
"\n",
" sch.pp.fit_all_sigmoids(adata, genes=adata.var['scHopfield_used'].values, spliced_key=SPLICED_KEY)\n",
" sch.pp.compute_sigmoid(adata, spliced_key=SPLICED_KEY)\n",
"\n",
" mse = adata.var.loc[adata.var['scHopfield_used'], 'sigmoid_mse']\n",
" print(f\"Sigmoid MSE: mean={mse.mean():.4f}, median={mse.median():.4f}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Scaffold: 90 TFs, 1857 targets, 75325 potential edges\n",
"TFs available for KO screening: 90\n"
]
}
],
"source": [
"# ── GRN scaffold from CellOracle mouse scATAC base GRN ───────────────────────\n",
"base_GRN = co.data.load_mouse_scATAC_atlas_base_GRN()\n",
"base_GRN.drop(['peak_id'], axis=1, inplace=True)\n",
"\n",
"gene_names = adata.var.index[adata.var['scHopfield_used'].values]\n",
"scaffold = pd.DataFrame(0, index=gene_names, columns=gene_names)\n",
"\n",
"tfs = list(set(base_GRN.columns.str.lower()) & set(scaffold.index.str.lower()))\n",
"targets = list(set(base_GRN['gene_short_name'].str.lower().values) & set(scaffold.columns.str.lower()))\n",
"index_map = {g.lower(): g for g in scaffold.index}\n",
"col_map = {g.lower(): g for g in scaffold.columns}\n",
"\n",
"for tf in tfs:\n",
" tf_col = [c for c in base_GRN.columns if c.lower() == tf][0]\n",
" for tgt in base_GRN[base_GRN[tf_col] == 1]['gene_short_name']:\n",
" if tgt.lower() in col_map:\n",
" scaffold.loc[index_map[tf], col_map[tgt.lower()]] = 1\n",
"\n",
"# Collect TF names (genes that appear as TFs in the scaffold)\n",
"TF_NAMES = sorted([index_map[tf] for tf in tfs])\n",
"print(f\"Scaffold: {len(tfs)} TFs, {len(targets)} targets, {int(scaffold.sum().sum())} potential edges\")\n",
"print(f\"TFs available for KO screening: {len(TF_NAMES)}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"if not LOAD_DATA:\n",
" # ── GRN inference ─────────────────────────────────────────────────────────────\n",
" sch.inf.fit_interactions(\n",
" adata,\n",
" cluster_key=CLUSTER_KEY,\n",
" spliced_key=SPLICED_KEY,\n",
" velocity_key=VELOCITY_KEY,\n",
" degradation_key=DEGRADATION_KEY,\n",
" n_epochs=N_EPOCHS,\n",
" batch_size=BATCH_SIZE,\n",
" device=DEVICE,\n",
" refit_gamma=True,\n",
" w_scaffold=scaffold.values.T,\n",
" scaffold_regularization=SCAFFOLD_REG,\n",
" reconstruction_regularization=100,\n",
" bias_regularization=1,\n",
" only_TFs=True,\n",
" w_threshold=1e-12,\n",
" skip_all=True,\n",
" learning_rate=0.1,\n",
" use_plateau_scheduler=True,\n",
" plateau_patience=100,\n",
" plateau_factor=0.1,\n",
" balanced_sampling=True,\n",
" drop_last=True,\n",
" include_neighbors=True,\n",
" neighbor_fraction=0.2,\n",
" get_plots=False,\n",
" )\n",
"\n",
" clusters = sorted(adata.obs[CLUSTER_KEY].unique())\n",
" print(f\"GRN inference complete for {len(clusters)} clusters.\")\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"if not LOAD_DATA:\n",
" # # ── Save adata checkpoint ─────────────────────────────────────────────────────\n",
" # # Save fitted model parameters using scHopfield's IO (W, I, gamma, sigmoid params).\n",
" sch.tl.save_model(adata, MODEL_PATH, overwrite=True)\n",
"\n",
" # # Save full adata (energies, Jacobians, centrality, UMAP, etc.).\n",
" # # Temporarily remove PyTorch model objects from uns — h5ad cannot serialise them.\n",
" # _models_bak = adata.uns.get('scHopfield', {}).pop('models', None)\n",
" # adata.write_h5ad(ADATA_PATH)\n",
" # if _models_bak is not None:\n",
" # adata.uns['scHopfield']['models'] = _models_bak\n",
"\n",
" # with open(TF_NAMES_PATH, 'w') as f:\n",
" # json.dump(TF_NAMES, f)\n",
"\n",
" # pd.DataFrame(\n",
" # {k: {'R': v[0], 'G': v[1], 'B': v[2], 'A': v[3]} for k, v in colors.items()}\n",
" # ).T.to_csv(COLORS_PATH)\n",
"\n",
" print(f\"Saved model → {MODEL_PATH}\")\n",
" # print(f\"Saved adata → {ADATA_PATH}\")\n",
" # print(f\"Saved TF names → {TF_NAMES_PATH} ({len(TF_NAMES)} TFs)\")\n",
" # print(f\"Saved colors → {COLORS_PATH}\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" #### ↳ Resume: load adata checkpoint\n",
"\n",
"\n",
"\n",
" Run this cell **instead of** the pipeline cells above (load data → velocity → sigmoid → scaffold → inference → energies → centrality → colours) when the checkpoint already exists.\n",
"\n",
" Sections 1–5 can then be run normally."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Model loaded from '/Users/bernaljp/Documents/SCHData/checkpoints/06_lineage_drivers/model.h5sch' | clusters=['10GMP', '11DC', '12Baso', '13Baso', '14Mo', '15Mo', '16Neu', '17Neu', '18Eos', '19Lymph', '1Ery', '2Ery', '3Ery', '4Ery', '5Ery', '6Ery', '7MEP', '8Mk', '9GMP'] | genes=1999\n",
"Loaded adata : 2671 cells × 1999 genes\n"
]
}
],
"source": [
"if LOAD_DATA:\n",
" # oracle_demo = co.data.load_tutorial_oracle_object()\n",
" # adata = oracle_demo.adata.copy()\n",
" sch.tl.load_model(adata, MODEL_PATH, overwrite=True)\n",
" # clusters = sorted(adata.obs[CLUSTER_KEY].unique())\n",
"\n",
" # with open(TF_NAMES_PATH) as f:\n",
" # TF_NAMES = json.load(f)\n",
"\n",
" # _col_df = pd.read_csv(COLORS_PATH, index_col=0)\n",
" # colors = {k: np.array([row.R, row.G, row.B, row.A]) for k, row in _col_df.iterrows()}\n",
"\n",
" print(f\"Loaded adata : {adata.n_obs} cells × {adata.n_vars} genes\")\n",
" # print(f\"Loaded TF names : {len(TF_NAMES)} TFs\")\n",
" # print(f\"Loaded colors : {len(colors)} clusters\")\n",
" adata\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Energies computed.\n",
"Network centrality computed.\n",
"Energy–gene correlations computed.\n"
]
}
],
"source": [
"# ── Compute energies, Jacobians, centrality, energy–gene correlation ───────────\n",
"sch.tl.compute_energies(adata, cluster_key=CLUSTER_KEY, spliced_key=SPLICED_KEY)\n",
"print(\"Energies computed.\")\n",
"\n",
"sch.tl.compute_network_centrality(adata, cluster_key=CLUSTER_KEY, threshold_number=40000)\n",
"print(\"Network centrality computed.\")\n",
"\n",
"sch.tl.energy_gene_correlation(adata, spliced_key=SPLICED_KEY, cluster_key=CLUSTER_KEY)\n",
"print(\"Energy–gene correlations computed.\")\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Colour map ready.\n"
]
}
],
"source": [
"# ── Extract colour map for consistent plots ───────────────────────────────────\n",
"fig_tmp, ax_tmp = plt.subplots()\n",
"scv.pl.scatter(adata, color=CLUSTER_KEY, basis=BASIS, ax=ax_tmp, show=False)\n",
"colors = {}\n",
"for k in adata.obs[CLUSTER_KEY].unique():\n",
" idx = np.where(adata.obs[CLUSTER_KEY] == k)[0][0]\n",
" c = ax_tmp.get_children()[0]._facecolors[idx].copy()\n",
" c[3] = 1.0\n",
" colors[k] = c\n",
"plt.close(fig_tmp)\n",
"print(\"Colour map ready.\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" ## 1. Candidate TF Prioritization\n",
"\n",
"\n",
"\n",
" We rank TFs by three complementary signals combined into a composite rank-sum score via\n",
"\n",
" `sch.tl.score_driver_tfs`:\n",
"\n",
"\n",
"\n",
" | Signal | Source | Description |\n",
"\n",
" |--------|--------|-------------|\n",
"\n",
" | W row L2-norm | `adata.varp['W_{cluster}']` | Regulatory strength |\n",
"\n",
" | Out-degree centrality | `adata.var['degree_centrality_out_{cluster}']` | Network hub-ness |\n",
"\n",
" | Energy–gene correlation | `adata.var['correlation_total_{cluster}']` | Expression drives energy landscape |\n",
"\n",
"\n",
"\n",
" **Composite score**: sum of ranks across signals per lineage;\n",
"\n",
" `lineage_bias = score_A − score_B` (erythroid minus myeloid)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# ── Score all TFs with the package API ────────────────────────────────────────\n",
"tf_df_all = sch.tl.score_driver_tfs(\n",
" adata,\n",
" lineage_A_clusters=ERYTHROID,\n",
" lineage_B_clusters=MYELOID,\n",
" cluster_key=CLUSTER_KEY,\n",
")\n",
"tf_df = tf_df_all.loc[[g for g in tf_df_all.index if g in TF_NAMES]].copy()\n",
"tf_df['is_TF'] = True\n",
"\n",
"# Rename to match downstream variable names (A = erythroid, B = myeloid)\n",
"tf_df = tf_df.rename(columns={\n",
" 'score_A': 'total_score_ery', 'score_B': 'total_score_mye',\n",
" 'wnorm_A': 'wnorm_ery', 'wnorm_B': 'wnorm_mye',\n",
" 'deg_A': 'deg_ery', 'deg_B': 'deg_mye',\n",
" 'ecorr_A': 'ecorr_ery', 'ecorr_B': 'ecorr_mye',\n",
"})\n",
"\n",
"TOP_N_INTERSECTION = 50\n",
"N_CANDIDATES_PER_LINEAGE = 8\n",
"\n",
"# ── Intersection-based exclusive candidate selection ──────────────────────────\n",
"top_ery_set = (\n",
" set(tf_df.nlargest(TOP_N_INTERSECTION, 'wnorm_ery').index) &\n",
" set(tf_df.nlargest(TOP_N_INTERSECTION, 'deg_ery').index) &\n",
" set(tf_df.nlargest(TOP_N_INTERSECTION, 'ecorr_ery').index)\n",
")\n",
"top_mye_set = (\n",
" set(tf_df.nlargest(TOP_N_INTERSECTION, 'wnorm_mye').index) &\n",
" set(tf_df.nlargest(TOP_N_INTERSECTION, 'deg_mye').index) &\n",
" set(tf_df.nlargest(TOP_N_INTERSECTION, 'ecorr_mye').index)\n",
")\n",
"\n",
"shared_ery_df = tf_df.loc[list(top_ery_set)].copy()\n",
"shared_mye_df = tf_df.loc[list(top_mye_set)].copy()\n",
"\n",
"exclusive_ery = [g for g in shared_ery_df.sort_values('total_score_ery', ascending=False).index if g not in top_mye_set]\n",
"exclusive_mye = [g for g in shared_mye_df.sort_values('total_score_mye', ascending=False).index if g not in top_ery_set]\n",
"\n",
"final_candidates_ery = exclusive_ery[:N_CANDIDATES_PER_LINEAGE]\n",
"final_candidates_mye = exclusive_mye[:N_CANDIDATES_PER_LINEAGE]\n",
"\n",
"final_ery_df = tf_df.loc[final_candidates_ery]\n",
"final_mye_df = tf_df.loc[final_candidates_mye]\n",
"\n",
"CANDIDATES = final_candidates_ery + final_candidates_mye\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"--- ERYTHROID: Top 8 Exclusive Candidates ---\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" total_score_ery | \n",
" total_score_mye | \n",
" wnorm_ery | \n",
" deg_ery | \n",
" ecorr_ery | \n",
"
\n",
" \n",
" \n",
" \n",
" | Ybx1 | \n",
" 5389.0 | \n",
" 5875.0 | \n",
" 0.312438 | \n",
" 0.166238 | \n",
" 0.432097 | \n",
"
\n",
" \n",
" | Stat1 | \n",
" 5248.0 | \n",
" 3399.0 | \n",
" 0.333781 | \n",
" 0.292721 | \n",
" 0.345669 | \n",
"
\n",
" \n",
" | Nr3c1 | \n",
" 4959.0 | \n",
" 4251.0 | \n",
" 0.203422 | \n",
" 0.461890 | \n",
" 0.483341 | \n",
"
\n",
" \n",
" | Myc | \n",
" 4920.0 | \n",
" 4246.0 | \n",
" 0.371766 | \n",
" 0.483698 | \n",
" 0.214282 | \n",
"
\n",
" \n",
" | Nfkb1 | \n",
" 4744.0 | \n",
" 4223.0 | \n",
" 0.339866 | \n",
" 0.208780 | \n",
" 0.191014 | \n",
"
\n",
" \n",
" | Rreb1 | \n",
" 4506.0 | \n",
" 4301.0 | \n",
" 0.216993 | \n",
" 0.228228 | \n",
" 0.293958 | \n",
"
\n",
" \n",
" | Zbtb7a | \n",
" 4500.0 | \n",
" 3803.0 | \n",
" 0.301138 | \n",
" 0.456885 | \n",
" 0.153372 | \n",
"
\n",
" \n",
" | Zbtb1 | \n",
" 4155.0 | \n",
" 4629.0 | \n",
" 0.236348 | \n",
" 0.257829 | \n",
" 0.157899 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" total_score_mye | \n",
" total_score_ery | \n",
" wnorm_mye | \n",
" deg_mye | \n",
" ecorr_mye | \n",
"
\n",
" \n",
" \n",
" \n",
" | Mef2c | \n",
" 5171.0 | \n",
" 2868.0 | \n",
" 0.248167 | \n",
" 0.334885 | \n",
" 0.399624 | \n",
"
\n",
" \n",
" | Irf2 | \n",
" 5163.0 | \n",
" 3826.0 | \n",
" 0.391317 | \n",
" 0.218168 | \n",
" 0.240946 | \n",
"
\n",
" \n",
" | Nfatc3 | \n",
" 4793.0 | \n",
" 4725.0 | \n",
" 0.219280 | \n",
" 0.210210 | \n",
" 0.337441 | \n",
"
\n",
" \n",
" | Irf8 | \n",
" 4775.0 | \n",
" 3525.0 | \n",
" 0.355602 | \n",
" 0.301301 | \n",
" 0.160587 | \n",
"
\n",
" \n",
" | Stat3 | \n",
" 4675.0 | \n",
" 4210.0 | \n",
" 0.248397 | \n",
" 0.563013 | \n",
" 0.229884 | \n",
"
\n",
" \n",
" | Runx1 | \n",
" 4582.0 | \n",
" 2962.0 | \n",
" 0.233669 | \n",
" 0.198599 | \n",
" 0.246073 | \n",
"
\n",
" \n",
" | Rel | \n",
" 4415.0 | \n",
" 4339.0 | \n",
" 0.253470 | \n",
" 0.183383 | \n",
" 0.172445 | \n",
"
\n",
" \n",
" | Irf1 | \n",
" 4262.0 | \n",
" 2883.0 | \n",
" 0.243509 | \n",
" 0.412212 | \n",
" 0.147591 | \n",
"
\n",
" \n",
"
\n",
"