MeRN Quick Start Tutorial

MeRN Quick Start Tutorial#

Metabolic Representation Net (MeRN)#

Single-cell RNA-seq data provides rich information about gene expression across thousands of cells, but interpreting cellular metabolism from transcriptomic data alone is challenging. Many metabolic processes are coordinated across genes and reactions, making them difficult to capture with standard dimensionality reduction methods.

Metabolic Representation Net (MeRN) is a graph-guided variational autoencoder designed to infer metabolic state from single-cell transcriptomic data.

Standard VAEs such as scVI learn low-dimensional latent representations of cells and decode them to gene-level parameters using flexible neural networks. While these models capture complex transcriptional structure, the latent dimensions are not explicitly linked to metabolic processes.

MeRN introduces a metabolic prior by leveraging a metabolic reaction graph, where:

  • nodes represent metabolic reactions

  • edges connect reactions that share metabolites

A graph VAE learns embeddings of this reaction network, and the cell-level latent space is decomposed into two components:

• Metabolic variation: variation explained by metabolic reactions
• Background variation: non-metabolic transcriptional variation

This structure links latent dimensions to reaction activity and enzyme-coding genes, enabling interpretable metabolic representations.

In this tutorial, we apply MeRN to a mouse embryogenesis single-cell RNA-seq dataset. We will:

  1. Load the preprocessed wild-type embryogenesis dataset.

  2. Fit or load a pretrained MeRN model.

  3. Extract metabolic and background cell embeddings.

  4. Analyze the metabolic latent space with UMAP and clustering.

  5. Decode reaction activity for every cell.

  6. Aggregate reaction activity into KEGG pathway scores.

  7. Build data-driven pathways (DDPs), which group graph-connected reactions with coordinated activity across cells.

  8. Compare reaction, KEGG pathway, and DDP activity across metabolic clusters.

  9. Visualize glycolysis/gluconeogenesis on KEGG pathway graphs, colored by DDP membership and by differential reaction activity.

  10. Expand a KEGG pathway by plotting all reactions in DDPs connected to glycolysis.

  11. Compare DDP structure between wild-type and Folr1 knockout embryos in a matched metabolic state. Steps 1-10 use the wild type alone; the knockout enters only in this step, as its own cells (embedded with the wild-type model to map them onto the wild-type metabolic states) and a precomputed table of its decoded reaction activity.

The tutorial emphasizes DDPs as a data-driven interpretation layer between individual reactions and broad canonical pathways. Reactions are highly granular, while KEGG pathways are predefined and often broad. DDPs provide coordinated metabolic reaction programs learned from the data and constrained by the metabolic graph.

MeRN Architecture

Files and Folder Structure#

tutorials/
├── MeRN_tutorial.ipynb
├── data/
│   ├── main_wt_anndata_pp_subsample.h5ad.gz     wild-type cells: log1p `X`, raw `counts_RNA`
│   ├── folr1_anndata_subsample.h5ad.gz          Folr1 knockout cells, same gene space
│   ├── folr1_average_enzyme_activity.csv.gz     knockout reaction activity
│   └── wt_ref_color_mapping.csv                 cell-type colors used in the paper's figures
├── figures/                                     figures written by the notebook
└── model/
    └── mern_rep_0/ … mern_rep_4/
        ├── model.pt.gz                          pretrained replicate checkpoints
        └── metabolic_clustering_0.4.csv         paper's metabolic clusters (mern_rep_0 only)

from scvi.train._callbacks import SaveCheckpoint
from mern import MERN
import mern.support as mern_support
import anndata as ad
import pandas as pd
import scanpy as sc
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import itertools
import os
import plotly.colors as pc
from contextlib import contextmanager
from IPython.display import display

train_new_model = False
clustering_res = 0.4
chosen_rep = 0
n_reps = 5

base_path_data = 'data'
base_path_model = 'model'
base_path_figures = 'figures'

wt_adata_file = 'main_wt_anndata_pp_subsample.h5ad.gz'
folr1_adata_file = 'folr1_anndata_subsample.h5ad.gz'
folr1_activity_file = 'folr1_average_enzyme_activity.csv.gz'
color_mapping_file = 'wt_ref_color_mapping.csv'

suffix = f'_rep{chosen_rep}'
metabolic_key = f"X_mern_metabolic{suffix}"
background_key = f"X_mern_background{suffix}"
metabolic_umap_key = f"X_umap_metabolic{suffix}"
background_umap_key = f"X_umap_background{suffix}"
graph_key = f"graph{suffix}"

Unpacking the Checkpoints#

The replicate checkpoints are gzipped and must be decompressed in place the first time the tutorial is run.

import gzip
import shutil
from pathlib import Path

for archive in sorted(Path(base_path_model).glob('mern_rep_*/model.pt.gz')):
    checkpoint = archive.with_suffix('')
    if checkpoint.exists():
        continue
    print(f'unpacking {archive}')
    with gzip.open(archive) as src, open(checkpoint, 'wb') as dst:
        shutil.copyfileobj(src, dst)

print('checkpoints ready:', sorted(p.parent.name for p in Path(base_path_model).glob('mern_rep_*/model.pt')))

checkpoints ready: ['mern_rep_0', 'mern_rep_1', 'mern_rep_2', 'mern_rep_3', 'mern_rep_4']

Figure Colors#

A small color registry keeps the figures in this tutorial consistent.

TEXT_COLOR = '#222222'
STRUCTURE_COLOR = '#666666'
BACKGROUND_COLOR = '#A7ADB8'
MISSING_COLOR = '#CFCFCF'
OBSERVED_COLOR = '#0072B2'
HIGHLIGHT_COLOR = '#D55E00'

METHOD_PALETTE = {
    'mern': '#0072B2',
    'rna': '#E69F00',
    'scrna': '#009E73',
    'prot': '#222222',
    'kegg_metroscreen': '#CC79A7',
    'compass': '#56B4E9',
}

GENOTYPE_PALETTE = {
    'WT': '#7A7A7A',
    'Folr1': '#8C6BB1',
    'Folr1 KO': '#8C6BB1',
}

# Cell-type colors as used in the manuscript figures
CELLTYPE_COLORS = pd.read_csv(f'{base_path_data}/{color_mapping_file}')
CELLTYPE_PALETTE = dict(zip(CELLTYPE_COLORS['cell_type'], CELLTYPE_COLORS['colors']))

GERMLAYER_PALETTE = {
    'Epiblast': '#000000',
    'PrimitiveStreak': '#808080',
    'GermCell': '#DB083D',
    'Endoderm': '#F5F200',
    'XMesoderm': '#00C4ED',
    'AxialMesoderm': '#0000FC',
    'Mesoderm': '#910299',
    'NeuralEctoderm': '#008C2A',
    'NonNeuralEctoderm': '#334513',
}

PHASE_PALETTE = {
    'G0/G1': '#4C78A8',
    'S': '#F2CF5B',
    'G2': '#72B7B2',
    'M': '#B279A2',
}

HEATMAP_CMAP = mpl.colormaps['viridis']

def reset_plotting_defaults():
    """Editable-text vector output and a consistent typeface."""
    mpl.rcdefaults()
    sns.reset_defaults()
    mpl.rcParams.update({
        'font.family': 'sans-serif',
        'font.sans-serif': ['Arial'],
        'axes.spines.top': False,
        'axes.spines.right': False,
        'legend.frameon': False,
        'figure.dpi': 150,
        'pdf.fonttype': 42,
        'ps.fonttype': 42,
        'svg.fonttype': 'none',
    })

reset_plotting_defaults()

def styled_theme(**kwargs):
    """`sns.set_theme` that keeps this tutorial's typeface and editable-text vector output."""
    kwargs.setdefault('font', 'Arial')
    rc = {
        'pdf.fonttype': 42,
        'ps.fonttype': 42,
        'svg.fonttype': 'none',
        'legend.frameon': False,
    }
    rc.update(kwargs.pop('rc', {}))
    sns.set_theme(rc=rc, **kwargs)

def top_features_per_cluster(z_scored, n_per_cluster):
    """The union of each cluster's most distinctive features, in nominating-cluster order.

    Ranking by variance across clusters says nothing about *which* cluster a feature belongs to, so
    a cluster whose features are all moderately elevated can go unrepresented. Taking each
    cluster's highest z-scores guarantees every cluster is characterised, and groups the rows into
    blocks that can be read against the columns.
    """
    picked = []
    for cluster in z_scored.columns:
        picked += z_scored[cluster].sort_values(ascending=False).head(n_per_cluster).index.tolist()
    return list(dict.fromkeys(picked))

def ddp_palette_sequence(labels):
    """Colors for a set of DDP labels, taken in order from Plotly's Dark24.

    DDP identifiers come from the clustering order, so they are specific to the cells a DDP set was
    built on. Colors are assigned by position rather than by identifier for that reason.
    """
    categories = list(pd.Series(labels).dropna().unique())
    palette = list(pc.qualitative.Dark24)
    return [palette[i % len(palette)] for i in range(len(categories))]

def ddp_color_map_for(labels):
    """DDP label -> color, for the `ddp_color_map` argument of the pathway plots."""
    categories = list(pd.Series(labels).dropna().unique())
    return dict(zip(categories, ddp_palette_sequence(categories)))

@contextmanager
def ddp_pathway_palette(labels):
    """The pathway plots color discrete edge labels from Plotly's Set1, in order of first
    appearance. Swap the DDP registry in for the duration of one plot."""
    original = pc.qualitative.Set1
    pc.qualitative.Set1 = ddp_palette_sequence(labels)
    try:
        yield
    finally:
        pc.qualitative.Set1 = original

Loading the WT Dataset#

The wild-type dataset here is a subsample of the original full embryogenesis data main_wt_anndata_pp.h5ad.gz. The pretrained models are the published ones.

Cells outside the embryo proper (the Blood&Vas, XEctoderm, and XEndoderm germ layers) and genes detected in fewer than 100 cells were already removed from the full object, so no further preprocessing is needed here.

The feature space is the object’s own highly_variable_metabolic / highly_variable_background flags, sorted (800 metabolic & 5,000 background genes). Commented code shows how to calculate highly variable genes for your own dataset.

X is already log1p-normalized; layers['counts_RNA'] holds the raw counts the model reads.

rna_full = ad.read_h5ad(f'{base_path_data}/{wt_adata_file}')

dataset = mern_support.KeggKGMLMetabolicDataset(species='mouse', capitalize=False)
dataset.add_module_info(rna_full)

# calculate highly variable genes yourself
"""temp_rna = rna_full[:, rna_full.var['Metabolic Gene'] == 'Metabolic']
sc.pp.highly_variable_genes(temp_rna, n_top_genes=800, flavor="seurat_v3", layer="counts_RNA")
sc.pl.highly_variable_genes(temp_rna, log=True)
met_hvgs = list(temp_rna.var_names[temp_rna.var["highly_variable"]])
temp_rna = rna_full[:, rna_full.var['Metabolic Gene'] == 'Non-metabolic']
sc.pp.highly_variable_genes(temp_rna, n_top_genes=5000, flavor="seurat_v3", layer="counts_RNA")
sc.pl.highly_variable_genes(temp_rna, log=True)
back_hvgs = list(temp_rna.var_names[temp_rna.var["highly_variable"]])
rna_full.var['highly_variable_metabolic'] = [True if gene in met_hvgs else False for gene in rna_full.var.index]
rna_full.var['highly_variable_background'] = [True if gene in back_hvgs else False for gene in rna_full.var.index]"""

features_to_use = sorted(
    rna_full.var_names[
        rna_full.var['highly_variable_metabolic'] | rna_full.var['highly_variable_background']
    ]
)
rna = rna_full[:, features_to_use].copy()

met_g = dataset.metabolic_topology(rna, self_loops=False)
dataset.add_rxn_module_info(rna, met_g)
dataset.add_rxn_module_info(rna_full, met_g)
rxn_genes = dataset.get_rxn_genes_all()

rna_full.layers['log1p_norm'] = rna_full.X
rna.layers['log1p_norm'] = rna.X

rxn_info = rna.uns["Reaction Info"]
rxn_ids = rxn_info.index.to_list()

print(rna)
print(rna.var['Metabolic Gene'].value_counts())

AnnData object with n_obs × n_vars = 6000 × 5800
    obs: 'nCount_RNA', 'nFeature_RNA', 'log_count', 'percent.mt', 'sex', 'embryo_id', 'sample_group', 'timepoint', 'cell_type', 'germlayer'
    var: 'names', 'n_cells', 'Gene Reactions', 'Metabolic Gene', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm', 'highly_variable_metabolic', 'highly_variable_background'
    uns: 'hvg', 'Reaction Info'
    obsm: 'X_pca', 'X_umap'
    layers: 'counts_RNA', 'log1p_norm'
Metabolic Gene
Non-metabolic    5000
Metabolic         800
Name: count, dtype: int64

for adata_obj in [rna, rna_full]:
    for col in ('cell_type', 'germlayer'):
        if col in adata_obj.obs:
            adata_obj.obs[col] = adata_obj.obs[col].astype('category')

    if 'germlayer' in adata_obj.obs:
        adata_obj.uns['germlayer_colors'] = [
            GERMLAYER_PALETTE.get(cat, MISSING_COLOR)
            for cat in adata_obj.obs['germlayer'].cat.categories
        ]

    if 'cell_type' in adata_obj.obs:
        adata_obj.uns['cell_type_colors'] = [
            CELLTYPE_PALETTE.get(cat, MISSING_COLOR)
            for cat in adata_obj.obs['cell_type'].cat.categories
        ]

Loading or Training a MeRN Model#

You can either train a new MeRN model or load the pretrained replicates. With train_new_model = True, the cell below sets the loss weights and KL schedule, checkpoints on the validation loss, trains with early stopping, and reloads the best checkpoint. Otherwise (the default), it loads the five pretrained replicate checkpoints: replicate chosen_rep supplies the embeddings, and all five are averaged when decoding reaction activity below.

if train_new_model:
    data_dir = 'models'

    plan_kwargs={
        "max_kl_weight": 0.001,
        "graph_kl_weight": 0.1,
        "data_elbo_weight": 1.0,
        "graph_elbo_weight": 0.2,
        "rxn_genes_weight": 0.5,
        "background_to_metabolic_weight": 30000,
        "lr": 0.002,
        "n_steps_kl_warmup": 0,
        "n_epochs_kl_warmup": None
    }

    checkpointing = SaveCheckpoint(
        dirpath=f'{data_dir}/scvi_log_3', monitor="validation_loss", load_best_on_end=False, check_nan_gradients=False
    )
    MERN.setup_anndata(rna, layer='counts_RNA', batch_key=None)
    mern = MERN(rna, met_g, rxn_genes, n_hidden=256, n_layers=2, n_metabolic_dim=25, n_background_dim=15, encode_covariates=False, strict_met_back_separation=True, rxn_genes_bias=True)
    mern.train(accelerator="auto", max_epochs=400, early_stopping=True, early_stopping_monitor="validation_loss", plan_kwargs=plan_kwargs, callbacks = [checkpointing])

    best_mern = MERN.load(checkpointing.best_model_path, rna, accelerator='auto')
    mern = best_mern
    models = [best_mern]
else:
    models_by_rep = {}
    for rep in range(n_reps):
        try:
            models_by_rep[rep] = MERN.load(
                f'{base_path_model}/mern_rep_{rep}', rna, accelerator='cpu'
            )
        except Exception as error:
            print(f'replicate {rep} unavailable ({type(error).__name__})')

    if chosen_rep not in models_by_rep:
        raise RuntimeError(f'replicate {chosen_rep} is required for the embeddings below')

    models = list(models_by_rep.values())
    mern = models_by_rep[chosen_rep]
    print(f'{len(models)} of {n_reps} replicates loaded')

INFO     File model/mern_rep_0/model.pt already downloaded                                                         
INFO     File model/mern_rep_1/model.pt already downloaded                                                         
INFO     File model/mern_rep_2/model.pt already downloaded                                                         
INFO     File model/mern_rep_3/model.pt already downloaded                                                         
INFO     File model/mern_rep_4/model.pt already downloaded                                                         
5 of 5 replicates loaded

Extracting the Embeddings#

The model maps every cell into two latent spaces (a metabolic one tied to the reaction graph, and a background one for everything else) and every reaction into a graph embedding. get_latent_representation returns all three. It reports the mean of each latent distribution, so repeated calls give identical answers.

metabolic_latent, background_latent, graph_embedding = mern.get_latent_representation(rna)

rna.obsm[metabolic_key] = metabolic_latent
rna.obsm[background_key] = background_latent
rna.uns[graph_key] = pd.DataFrame(graph_embedding, index=mern.vertex_names_ordered)

print(f'metabolic latent : {metabolic_latent.shape[0]:,} cells x {metabolic_latent.shape[1]} dims')
print(f'background latent: {background_latent.shape[0]:,} cells x {background_latent.shape[1]} dims')
print(f'graph embedding  : {graph_embedding.shape[0]:,} reactions x {graph_embedding.shape[1]} dims')

metabolic latent : 6,000 cells x 25 dims
background latent: 6,000 cells x 15 dims
graph embedding  : 1,213 reactions x 25 dims

Metabolic and Background Embeddings#

MeRN learns two latent embeddings for each cell:

Background embedding#

The background embedding captures non-metabolic transcriptional variation. This can include broad cell identity, developmental programs, technical structure, and other gene expression patterns that are not specifically explained through the metabolic reaction graph.

Metabolic embedding#

The metabolic embedding captures variation that is routed through the metabolic reaction graph. Cells that are transcriptionally similar can still separate in this space if they differ in inferred reaction activity or coordinated metabolic programs.

This embedding is the main space we use to identify metabolic cell states.

By comparing the metabolic and background embeddings, we can ask whether differences between cells are better explained by:

  • general transcriptional variation

  • metabolic reaction-level variation

In the next cells, we compute neighborhoods, UMAPs, and clusters separately for the metabolic and background embeddings.

Embedding Clustering#

sc.pp.neighbors(rna, use_rep=metabolic_key, key_added='met_neighbors')
sc.tl.umap(rna, neighbors_key='met_neighbors')
sc.tl.leiden(rna, neighbors_key='met_neighbors', key_added='met_leiden', resolution=clustering_res)
sc.pl.umap(rna, color=['met_leiden', 'cell_type'], neighbors_key='met_neighbors')
../_images/7a4de9d0c03872ca7128af5fab3e2cd73f2c7afd56d2e596a3708e9f2a4424af.png

Because this tutorial uses a subsample of the full dataset, the neighbor graph, UMAP, and Leiden clustering are recomputed on fewer cells. Thus, the UMAP layout differs from the paper’s figure, and Leiden partitions the cells slightly differently. We load in the metabolic clusters from the paper below to offer a side by side comparison, but the downstream analyses in this tutorial use the newly computed met_leiden clusters.

paper_clust_path = './model/mern_rep_0/metabolic_clustering_0.4.csv'
clust = pd.read_csv(paper_clust_path, index_col=0)['Metabolic Clustering_0_0.4']

missing = rna.obs_names.difference(clust.index)
assert len(missing) == 0, f'{len(missing)} cells in rna not found in the clustering CSV, e.g. {list(missing[:5])}'

rna.obs['paper_met_clust'] = clust.reindex(rna.obs_names).astype(str).astype('category')
sc.pl.umap(rna, color=['paper_met_clust', 'met_leiden', 'cell_type'], neighbors_key='met_neighbors')
../_images/b3186bec4348f51cfd267939d8f2757e2254fa24766b351c2cd24b312b8b549d.png 
sc.pp.neighbors(rna, use_rep=background_key, key_added='back_neighbors')
sc.tl.umap(rna, neighbors_key='back_neighbors')
sc.tl.leiden(rna, neighbors_key='back_neighbors', key_added='back_leiden', resolution=clustering_res)
sc.pl.umap(rna, color=['back_leiden', 'germlayer'], neighbors_key='back_neighbors')
../_images/726081ca911d32e6484464552b09f299c9e7f9f7f1a23c1bd670c7b893288838.png

Metabolic Embedding Analysis#

We cluster cells in the MeRN metabolic embedding to identify metabolic states. To relate these states to embryonic development, we first plot the germ-layer composition of each metabolic cluster, then the proportion of cells in each cluster across timepoints.

cluster_key = "met_leiden"
group_key = "germlayer"

rna.obs[cluster_key] = rna.obs[cluster_key].astype(str).astype("category")

counts = pd.crosstab(
    rna.obs[cluster_key],
    rna.obs[group_key],
)

props = counts.div(counts.sum(axis=1), axis=0)

germlayer_order = list(rna.obs[group_key].cat.categories)
props = props.reindex(columns=germlayer_order, fill_value=0)

colors = [
    GERMLAYER_PALETTE.get(germlayer, MISSING_COLOR)
    for germlayer in germlayer_order
]

ax = props.plot(
    kind="bar",
    stacked=True,
    color=colors,
    figsize=(8, 5),
    width=0.85,
    edgecolor="none",
)

ax.set_xlabel("Metabolic Cluster")
ax.set_ylabel("Proportion of Cells")
ax.set_title("Germ Layer Composition by Metabolic Cluster")
ax.set_ylim(0, 1)
ax.tick_params(axis="x", rotation=0)
ax.legend(
    title="Germ Layer",
    bbox_to_anchor=(1.02, 1),
    loc="upper left",
    frameon=False,
)

plt.tight_layout()
plt.show()
../_images/7748fd4bb6a7a5a47cb2ec3ea8afa7ad79baf36d3e36e203e60d4772122b1d8a.png 
time_key = "timepoint"
cluster_key = "met_leiden"

rna.obs[time_key] = rna.obs[time_key].astype(str).astype("category")
rna.obs[cluster_key] = rna.obs[cluster_key].astype(str).astype("category")

counts = pd.crosstab(
    rna.obs[time_key],
    rna.obs[cluster_key],
)

props = counts.div(counts.sum(axis=1), axis=0)

cluster_order = sorted(
    rna.obs[cluster_key].cat.categories,
    key=lambda x: int(x) if str(x).isdigit() else str(x),
)

props = props.reindex(columns=cluster_order, fill_value=0)

cluster_palette = dict(
    zip(
        cluster_order,
        sns.color_palette("tab10", n_colors=len(cluster_order)).as_hex(),
    )
)

colors = [cluster_palette[c] for c in cluster_order]

ax = props.plot(
    kind="bar",
    stacked=True,
    color=colors,
    figsize=(8, 5),
    width=0.85,
    edgecolor="none",
)

ax.set_xlabel("Timepoint")
ax.set_ylabel("Proportion of Cells")
ax.set_title("Metabolic Cluster Composition Over Time")
ax.set_ylim(0, 1)
ax.tick_params(axis="x", rotation=0)
ax.legend(
    title="Metabolic Cluster",
    bbox_to_anchor=(1.02, 1),
    loc="upper left",
    frameon=False,
)

plt.tight_layout()
plt.show()
../_images/0f21edc4fffa0eb478a091c21d6b7426446692c8cf4246a3e15aed7727b28da6.png

Reaction Activity, KEGG Pathway, and DDP Calculations#

MeRN decodes each cell into reaction activity: how active every metabolic reaction is inferred to be in that cell. From this reaction activity matrix we build two coarser summaries — canonical KEGG pathway scores, and data-driven pathways (DDPs) learned from the data itself.

That gives three activity matrices, at increasing levels of aggregation, which we use throughout the rest of the tutorial:

  1. Reaction activity: cells × reactions

  2. KEGG pathway activity: cells × canonical KEGG pathways

  3. DDP activity: cells × data-driven pathways

Every later section presents these three in the same order, from the finest to the coarsest view.

A note on naming: the package calls this quantity enzyme activity in its function and variable names (mern_support.average_enzyme_activity, enzyme_acts_mean, and the mern average-enzyme-activity CLI command). It is the same matrix. This tutorial says reaction activity in the text and figures, because the values are indexed by reaction rather than by gene.

DDPs

Decoding Reaction Activity#

MeRN decodes each cell into reaction-level activity: a cells x reactions matrix giving the activity inferred for each reaction in each cell.

Decoding is stochastic, and separately trained replicates differ from one another, so the activity used for analysis is averaged rather than taken from a single decode. mern_support.average_enzyme_activity handles both: it decodes n_samples times per model and averages across every replicate it is given.

To demonstrate, we average 8 decodes from each of the 5 included replicates below.

enzyme_acts_mean = mern_support.average_enzyme_activity(models, rna, n_samples=8, show_progress=False)

enzyme_acts_mean

Seed set to 8
Seed set to 9
Seed set to 10
Seed set to 11
Seed set to 12
Seed set to 13
Seed set to 14
Seed set to 15
Seed set to 1008
Seed set to 1009
Seed set to 1010
Seed set to 1011
Seed set to 1012
Seed set to 1013
Seed set to 1014
Seed set to 1015
Seed set to 2008
Seed set to 2009
Seed set to 2010
Seed set to 2011
Seed set to 2012
Seed set to 2013
Seed set to 2014
Seed set to 2015
Seed set to 3008
Seed set to 3009
Seed set to 3010
Seed set to 3011
Seed set to 3012
Seed set to 3013
Seed set to 3014
Seed set to 3015
Seed set to 4008
Seed set to 4009
Seed set to 4010
Seed set to 4011
Seed set to 4012
Seed set to 4013
Seed set to 4014
Seed set to 4015
rn:R00351 rn:R00352 rn:R00209 rn:R00212 rn:R01196 rn:R10866 rn:R00341 rn:R00345 rn:R00346 rn:R00431 rn:R00726 rn:R00217 rn:R00344 rn:R12212 rn:R00705 rn:R00706 rn:R00233 rn:R00353 rn:R00742 rn:R00357 rn:R03858 rn:R04742 rn:R03778 ... rn:R13391 rn:R05923 rn:R05924 rn:R13395 rn:R13402 rn:R13408 rn:R08781 rn:R02161 rn:R13224 rn:R00081
WT65_1_CTCCCAATCATCAGTG-1 0.501264 1.423109 2.900505 3.117703 0.000221 0.068973 2.125863 0.016358 0.000005 0.023986 ... 0.319463 0.142891 0.198835 1.290856 0.185172 0.681304 0.739687 0.836421 0.744113 0.715358
WT65_1_TCTAACTGTCATCACA-1 1.279655 3.017212 3.519320 4.090299 0.337239 0.263835 1.991866 0.169117 0.341080 0.383189 ... 0.316503 0.226324 0.325084 1.498919 0.420310 0.674704 0.971236 0.957889 0.720279 0.787069
WT65_1_AGCATCACACTGCACG-1 0.725717 2.303315 3.178614 3.547062 0.339189 0.098212 1.816407 0.094003 0.133307 0.164853 ... 0.313879 0.073306 0.220451 1.530750 0.208325 0.738438 0.780331 0.806572 0.639561 0.935785
WT65_1_AAGAACAAGAATCGCG-1 1.099178 2.751644 4.373340 4.908928 0.161633 0.194620 3.403841 0.030333 0.129293 0.221094 ... 0.345708 0.251665 0.309978 1.668479 0.138314 0.823985 1.123434 1.070077 0.885505 0.925709
WT65_1_AGACAAATCCAACCGG-1 2.997977 5.055899 4.469263 5.092192 1.186582 0.907151 2.681515 0.951997 1.271756 1.186453 ... 0.424448 0.188741 0.228933 1.538487 0.193824 0.888587 0.930118 0.985039 0.904613 1.042986
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
WT90_4_GTAGTACTCCGCGATG-1 3.265272 3.250404 3.044608 2.534521 1.608961 1.359930 1.974718 1.303628 1.182495 1.662586 ... 0.904989 1.096877 1.089835 0.571538 0.570002 0.370623 0.285955 0.521686 0.706695 0.721114
WT90_4_GCCGATGCACGGGCTT-1 0.677042 0.701352 0.686832 0.687992 0.317901 0.972473 0.179103 0.698784 0.644020 0.988871 ... 0.899537 0.840459 0.962301 0.739031 0.671718 1.796167 2.410409 2.895690 2.117692 2.157468
WT90_4_TCGTGCTAGAGCCGAT-1 1.484363 1.678304 2.925100 2.490174 0.287045 0.450374 1.648407 0.403327 0.551135 0.402139 ... 0.421659 1.402514 1.293007 0.570563 1.221698 2.419944 3.681092 4.067152 2.563948 3.034113
WT90_4_CAGATACCACCCTTGT-1 4.127961 3.907559 1.786773 1.759523 3.567367 3.800569 2.153339 3.464602 3.574801 3.677756 ... 0.433399 0.189135 0.400398 0.490558 0.708944 1.966772 2.870537 2.760200 2.104390 2.124063
WT90_4_GTCTCACGTTGCTCAA-1 6.363608 5.642426 2.614729 2.475686 3.928648 4.964446 1.524732 4.931694 4.693810 4.977001 ... 0.633721 0.667833 0.872392 0.573266 0.315242 1.506813 1.945249 2.141428 1.801134 1.973704

6000 rows × 1213 columns

Calculating KEGG Pathway Scores#

The first summary of reaction activity is the canonical one: KEGG’s own pathway definitions.

The input reaction activity matrix contains cells as rows and metabolic reactions as columns. Using the reaction metadata in rxn_info, each reaction is mapped to one or more KEGG pathways. For each cell and each KEGG pathway, we calculate the average reaction activity over the reactions assigned to that pathway.

This produces three outputs:

  • pathway_scores: a cells × KEGG pathways matrix of pathway activity scores

  • pathway_rxns: a dictionary mapping each KEGG pathway ID to its member reactions

  • pathway_names: a dictionary mapping KEGG pathway IDs to readable pathway names

These pathway scores give a canonical pathway-level view of metabolic activity. Because the pathway boundaries are drawn in advance, they are broad and identical for every dataset, which is what motivates the data-driven alternative in the next section.

pathway_scores, pathway_rxns, pathway_names = mern_support.calculate_pathway_scores(
    enzyme_acts=enzyme_acts_mean,
    reaction_info=rxn_info
)

print("\nKEGG pathway activity matrix")
print("  cells x pathways:", pathway_scores.shape)

pathway_sizes = pd.Series(
    {pathway: len(rxns) for pathway, rxns in pathway_rxns.items() if pathway in pathway_scores.columns},
    name="n_reactions",
)

display(pathway_sizes.describe())
display(pathway_scores.head())

KEGG pathway activity matrix
  cells x pathways: (6000, 122)

count     122.000000
mean       30.229508
std       112.081197
min         1.000000
25%         4.000000
50%        14.000000
75%        27.000000
max      1213.000000
Name: n_reactions, dtype: float64
rn00020 rn00630 rn01100 rn01110 rn01120 rn01200 rn01210 rn01230 rn00720 rn00620 ... rn00604 rn00625 rn00930 rn00740 rn00750 rn00830 rn00232 rn00660 rn00524 rn00531
WT65_1_CTCCCAATCATCAGTG-1 0.992666 1.952201 0.925399 1.300436 1.413782 1.972407 1.087445 2.615239 0.920727 1.228447 ... 0.181140 0.894273 0.104186 1.731454 0.289098 0.219887 0.121115 0.442724 0.178332 0.295616
WT65_1_TCTAACTGTCATCACA-1 1.324880 2.187607 0.985935 1.319960 1.613052 2.307511 1.240203 2.776433 1.310102 1.641467 ... 0.196100 0.970507 0.136696 1.727014 0.271104 0.145457 0.071652 0.438612 0.398476 0.348757
WT65_1_AGCATCACACTGCACG-1 1.056939 1.905424 1.010693 1.464017 1.579590 2.194766 1.146654 2.822768 0.987082 1.446957 ... 0.189852 0.965385 0.134730 1.632272 0.462235 0.301350 0.139733 0.402972 0.317831 0.241790
WT65_1_AAGAACAAGAATCGCG-1 1.567621 2.118088 1.037528 1.471715 1.756105 2.407500 1.442160 3.107196 1.402735 1.895322 ... 0.301406 0.915207 0.169464 1.479127 0.418618 0.211665 0.261326 0.475032 0.615253 0.274273
WT65_1_AGACAAATCCAACCGG-1 1.804183 2.496473 1.166687 1.685167 2.035642 2.912792 1.603595 3.328616 1.812392 2.318196 ... 0.151014 1.013375 0.180698 1.470742 0.317617 0.357161 0.215039 0.559147 0.331385 0.308630

5 rows × 122 columns

Calculating Data-Driven Pathways (DDPs)#

KEGG pathways are defined in advance and are the same for every dataset. DDPs are the data-driven counterpart: groups of reactions that are both connected in the metabolic reaction graph and coordinated across the cells of this dataset.

To construct DDPs, we:

  1. Compute reaction-reaction correlations within the tutorial’s epiblast-like metabolic cluster, selected by its Epiblast + PrimitiveStreak fraction.

  2. Cluster graph-connected reactions using a minimum correlation threshold.

  3. Score each cell for each DDP by averaging reaction activity across member reactions.

The resulting ddp_scores matrix has all WT cells as rows and DDPs as columns. We retain this state’s WT correlations and linkage for the later weakest-link comparison; the WT DDPs are calculated only here.

ddp_min_corr = 0.7
ddp_min_size = 3

# Use the same tutorial metabolic state for DDP construction and the KO comparison.
wt_cluster_labels = rna.obs['met_leiden'].astype(str)
germ_fractions = pd.crosstab(
    wt_cluster_labels, rna.obs['germlayer'].astype(str), normalize='index'
)
epiblast_like = [c for c in ('Epiblast', 'PrimitiveStreak') if c in germ_fractions.columns]
epiblast_score = germ_fractions[epiblast_like].sum(axis=1).sort_values(ascending=False)
comparison_state = str(epiblast_score.index[0])
wt_state_cells = rna.obs_names[wt_cluster_labels == comparison_state]

print('epiblast + primitive streak fraction by metabolic state:')
print(epiblast_score.round(3).to_string())
print(f'WT DDPs use tutorial state {comparison_state}: {len(wt_state_cells)} cells')

wt_corr = enzyme_acts_mean.loc[wt_state_cells].corr(method='spearman')

rxn_to_ddp, ddp_rxns, ddp_clusters, wt_linkage = mern_support.calculate_ddps(
    G=met_g,
    corr_df=wt_corr,
    min_corr=ddp_min_corr,
    min_size=ddp_min_size,
)

print("\nDDPs")
print("  number of DDPs:", len(ddp_rxns))

ddp_sizes = pd.Series(
    {ddp: len(rxns) for ddp, rxns in ddp_rxns.items()},
    name="n_reactions",
)

display(ddp_sizes.describe())
display(ddp_sizes.sort_values(ascending=False).head(10))

ddp_scores, ddp_rxns = mern_support.calculate_ddp_scores(
    enzyme_acts=enzyme_acts_mean,
    rxn_to_ddp=rxn_to_ddp,
)

print("\nDDP activity matrix")
print("  cells x DDPs:", ddp_scores.shape)

epiblast + primitive streak fraction by metabolic state:
met_leiden
2    0.924
4    0.046
3    0.042
1    0.035
0    0.024
WT DDPs use tutorial state 2: 1040 cells

DDPs
  number of DDPs: 115

DDP activity matrix
  cells x DDPs: (6000, 115)

count    115.000000
mean       7.730435
std        6.966488
min        3.000000
25%        3.000000
50%        5.000000
75%        9.000000
max       36.000000
Name: n_reactions, dtype: float64

ddp_114    36
ddp_113    34
ddp_112    33
ddp_111    27
ddp_110    23
ddp_109    23
ddp_108    22
ddp_107    20
ddp_106    20
ddp_105    19
Name: n_reactions, dtype: int64

Metabolic Activity Summaries Across Clusters#

We next summarize MeRN-derived metabolic activity across the metabolic clusters identified from the metabolic embedding.

We compare three levels of metabolic organization, in the same order as the previous section:

  1. Reaction activity: individual MeRN-inferred reaction activities

  2. KEGG pathway activity: canonical pathway-level summaries

  3. DDP activity: data-driven metabolic programs built from coordinated graph-connected reactions

For each matrix, we z-score each feature across all cells and then average those z-scores within each metabolic cluster. This highlights features that are relatively higher or lower in specific metabolic states.

Reaction-Level Activity Across Metabolic Clusters#

We first visualize reaction-level activity. Each row is a metabolic reaction and each column is a metabolic cluster. Each reaction is z-scored across all cells, and the values shown are the mean of those cell-level z-scores within each cluster.

There are too many reactions to show at once, so each cluster contributes 15 reactions with the highest mean z-score in that cluster.

styled_theme(style="white", font_scale=1.0)

cluster_key = "met_leiden"
cluster_series = rna.obs[cluster_key].astype(str)

acts = enzyme_acts_mean.loc[:, enzyme_acts_mean.std(axis=0) > 0]

cell_z = acts.sub(acts.mean(axis=0), axis=1).div(acts.std(axis=0), axis=1)
z_scored = cell_z.groupby(cluster_series, observed=False).mean().T

n_rxns_per_cluster = 15
plot_data = z_scored.loc[top_features_per_cluster(z_scored, n_rxns_per_cluster)]
plot_data.index = plot_data.index.astype(str)
print(f'{len(plot_data)} reactions from {n_rxns_per_cluster} per cluster')


def rxn_display_name(rxn_id, max_chars=44):
    """A readable name for a reaction node.

    `Rxn Name` holds the KEGG enzyme name of every reaction in the node, joined by ';' where a
    node groups several; the first is enough to identify the chemistry on a figure.
    """
    raw = str(rxn_info['Rxn Name'].get(rxn_id, '')).split(';')[0].strip()
    if not raw:
        raw = rxn_id
    return raw if len(raw) <= max_chars else raw[:max_chars - 1] + '…'


# Labelling every row is unreadable, so call out a defined selection instead of none: the two
# most cluster-specific reactions per metabolic cluster, plus the most variable reactions overall.
n_per_cluster = 2
n_overall = 4
callout_rxns = []
for cluster in plot_data.columns:
    callout_rxns += (
        plot_data[cluster].sort_values(ascending=False).head(n_per_cluster).index.tolist()
    )
callout_rxns += plot_data.var(axis=1).sort_values(ascending=False).head(n_overall).index.tolist()
callout_rxns = list(dict.fromkeys(callout_rxns))
print(f'calling out {len(callout_rxns)} of the {len(plot_data)} reactions shown')

g = sns.clustermap(
    plot_data,
    cmap=HEATMAP_CMAP,
    center=0,
    row_cluster=False,
    col_cluster=False,
    xticklabels=True,
    yticklabels=False,
    figsize=(9, 10),
    dendrogram_ratio=(0.08, 0.02),
    cbar_pos=(0.90, 0.28, 0.022, 0.22),
)

g.ax_heatmap.set_xlabel("Metabolic Cluster", fontsize=11)

g.ax_heatmap.set_ylabel("")
g.ax_heatmap.tick_params(axis="x", rotation=0)

g.fig.suptitle("Cluster-Averaged Reaction Activity", fontsize=14, y=0.985)
g.fig.text(
    0.24, 0.955, f"top {n_rxns_per_cluster} reactions per cluster, mean cell z-score",
    ha="center", fontsize=9, color=STRUCTURE_COLOR,
)
g.fig.subplots_adjust(left=0.04, right=0.44, top=0.94, bottom=0.06)
g.cax.set_position([0.90, 0.28, 0.022, 0.22])

n_rows, n_cols = plot_data.shape
row_order = (
    list(g.dendrogram_row.reordered_ind)
    if g.dendrogram_row is not None
    else list(range(n_rows))
)
display_pos = {plot_data.index[original]: drawn for drawn, original in enumerate(row_order)}

points = sorted((display_pos[rxn] + 0.5, rxn) for rxn in callout_rxns)

min_gap = n_rows / (len(points) + 1) * 0.78
label_y = []
for y, _ in points:
    label_y.append(y if not label_y else max(y, label_y[-1] + min_gap))
overflow = label_y[-1] - (n_rows - 0.5)
if overflow > 0:
    label_y = [y - overflow for y in label_y]

for (row_y, rxn), text_y in zip(points, label_y):
    g.ax_heatmap.annotate(
        rxn_display_name(rxn),
        xy=(n_cols, row_y),
        xytext=(n_cols + 0.9, text_y),
        textcoords="data",
        fontsize=7.5,
        va="center",
        ha="left",
        color=TEXT_COLOR,
        arrowprops=dict(
            arrowstyle="-",
            color=STRUCTURE_COLOR,
            linewidth=0.6,
            shrinkA=0,
            shrinkB=2,
        ),
        annotation_clip=False,
    )

plt.show()

75 reactions from 15 per cluster
calling out 12 of the 75 reactions shown
../_images/79729d31c287205df634aad0383741f427cc6eb13843f8fccfc5b33bb59585f1.png

KEGG Pathway Activity Across Metabolic Clusters#

Next, we aggregate reaction activity into KEGG pathway scores. Each row is a KEGG pathway and each column is a metabolic cluster. Each pathway score is z-scored across all cells, and the values shown are the mean of those cell-level z-scores within each cluster.

Similar to the reaction-level plot above, each cluster contributes the 5 pathways it scores highest for.

styled_theme(style="white", font_scale=1.0)

cluster_key = "met_leiden"
cluster_series = rna.obs[cluster_key].astype(str)

exclude_pathways = {"rn01100", "rn01110", "rn01120", "rn01200"}

ps = pathway_scores.loc[:, ~pathway_scores.columns.isin(exclude_pathways)]
ps = ps.loc[:, ps.std(axis=0) > 0]

cell_z = ps.sub(ps.mean(axis=0), axis=1).div(ps.std(axis=0), axis=1)
z_scored = cell_z.groupby(cluster_series, observed=False).mean().T

n_pathways_per_cluster = 5
plot_data = z_scored.loc[top_features_per_cluster(z_scored, n_pathways_per_cluster)]
print(f'{len(plot_data)} pathways from {n_pathways_per_cluster} per cluster')

plot_data.index = [pathway_names.get(pid, pid) for pid in plot_data.index]

g = sns.clustermap(
    plot_data,
    cmap=HEATMAP_CMAP,
    center=0,
    row_cluster=False,
    col_cluster=False,
    xticklabels=True,
    yticklabels=True,
    figsize=(10, 9),
    dendrogram_ratio=(0.12, 0.05),
    cbar_pos=(0.97, 0.2, 0.03, 0.45),
)

plt.setp(g.ax_heatmap.get_xticklabels(), rotation=0, fontsize=10)
plt.setp(g.ax_heatmap.get_yticklabels(), rotation=0, fontsize=9)

g.ax_heatmap.set_xlabel("Metabolic Cluster", fontsize=12)
g.ax_heatmap.set_ylabel("KEGG Pathway", fontsize=12)
g.fig.suptitle(
    f"Cluster-Averaged KEGG Pathway Activity\n"
    f"top {n_pathways_per_cluster} pathways per cluster, mean cell z-score",
    fontsize=14, y=1.04,
)

plt.show()

25 pathways from 5 per cluster
../_images/2e9558693bf8d32f357fbe9fe5d987a6c6f889d1987be1dbc86625ed0ffd081b.png

DDP Activity Across Metabolic Clusters#

Finally, we summarize activity at the DDP level. Each row is a data-driven pathway and each column is a metabolic cluster. Each DDP score is z-scored across all cells, and the values shown are the mean of those cell-level z-scores within each cluster.

DDPs provide an intermediate view between individual reactions and broad KEGG pathways as they group reactions that are connected in the metabolic graph and coordinated across cells.

styled_theme(style="white", font_scale=1.0)

cluster_key = "met_leiden"
cluster_series = rna.obs[cluster_key].astype(str)

ds = ddp_scores.loc[:, ddp_scores.std(axis=0) > 0]

cell_z = ds.sub(ds.mean(axis=0), axis=1).div(ds.std(axis=0), axis=1)
z_scored = cell_z.groupby(cluster_series, observed=False).mean().T

n_ddps_per_cluster = 5
plot_data = z_scored.loc[top_features_per_cluster(z_scored, n_ddps_per_cluster)]
print(f'{len(plot_data)} DDPs from {n_ddps_per_cluster} per cluster')

plot_data.index = plot_data.index.astype(str)

g = sns.clustermap(
    plot_data,
    cmap=HEATMAP_CMAP,
    center=0,
    row_cluster=False,
    col_cluster=False,
    xticklabels=True,
    yticklabels=True,
    figsize=(6.5, 7),
    dendrogram_ratio=(0.16, 0.02)
)

g.ax_heatmap.set_xlabel("Metabolic Cluster", fontsize=11)
g.ax_heatmap.set_ylabel("DDP", fontsize=11)
g.ax_heatmap.tick_params(axis="x", rotation=0)
g.ax_heatmap.tick_params(axis="y", labelsize=8)

g.fig.suptitle(
    f"Cluster-Averaged DDP Activity\n"
    f"top {n_ddps_per_cluster} DDPs per cluster, mean cell z-score",
    fontsize=13, y=1.0,
)

g.fig.subplots_adjust(left=0.08, right=0.72, top=0.94, bottom=0.08)
g.cax.set_position([0.88, 0.30, 0.035, 0.25])

plt.show()

25 DDPs from 5 per cluster
../_images/0b72cb1d956f67cf746a85086bea6942c2ff181ea0a5fdc62126e14a69cf0a05.png

Visualizing Metabolic Programs on KEGG Pathway Graphs#

After summarizing reaction, KEGG pathway, and DDP activity across metabolic clusters, we can project these signals back onto metabolic pathway structure. This helps connect the matrix-level summaries to specific reactions and metabolites within known biochemical pathways.

Here we use glycolysis/gluconeogenesis as an example pathway.

Glycolysis Reactions Colored by DDP Membership#

In this plot, each edge represents a reaction in the KEGG glycolysis/gluconeogenesis pathway (rn00010). We color each reaction by the DDP it belongs to.

pathway_id_to_plot = "rn00010"

ddp_edge_labels = rxn_to_ddp.copy()

with ddp_pathway_palette(ddp_edge_labels):
    fig = mern_support.kegg_pathway_plot(
        pathway=pathway_id_to_plot,
        dataset=dataset,
        rna=rna,
        user_labels=ddp_edge_labels,
        edge_width=4,
        show_node_labels=True,
    )

fig.update_layout(
    title="Glycolysis / Gluconeogenesis: Reactions Colored by DDP Membership"
)

fig.write_image(f"{base_path_figures}/glycolysis_ddp_plot.png", width=1200, height=800, scale=2)
fig.show()

Downloading from KEGG for mouse: rn00010. This uses KEGG access provided for academic use by academic users; you are responsible for complying with KEGG's terms. the result will be cached in memory for this dataset instance and reused.

Glycolysis DDP plot

Glycolysis Reactions Colored by Differential Activity#

We next color the glycolysis pathway by reaction-level differential activity between two metabolic clusters. Leiden numbers its clusters in an order that depends on the cells present, so instead of naming two ids this cell picks the contrast the figure is meant to show: the metabolic cluster with the highest mean glycolysis activity against the one with the lowest, among clusters large enough for the test to mean anything.

For each reaction, we compare inferred reaction activity between those two clusters using a Wilcoxon rank-sum test. The edge color shows Cohen’s d effect size:

  • positive values indicate higher inferred reaction activity in the high-glycolysis cluster

  • negative values indicate higher inferred reaction activity in the low-glycolysis cluster

pathway_id_to_plot = "rn00010"

cluster_key = "met_leiden"
cluster_labels = rna.obs[cluster_key].astype(str)

min_cluster_cells = 100
assert pathway_id_to_plot in pathway_scores.columns, 'glycolysis is not in the pathway scores'

cluster_sizes = cluster_labels.value_counts()
eligible = cluster_sizes.index[cluster_sizes >= min_cluster_cells]
glycolysis_by_cluster = (
    pathway_scores[pathway_id_to_plot]
    .groupby(cluster_labels, observed=False)
    .mean()
    .reindex(eligible)
    .dropna()
)
up_cluster = str(glycolysis_by_cluster.idxmax())
down_cluster = str(glycolysis_by_cluster.idxmin())
print(f'cluster {up_cluster} (highest mean glycolysis activity) vs '
      f'cluster {down_cluster} (lowest), of {len(eligible)} clusters with '
      f'at least {min_cluster_cells} cells')

up_cells = rna.obs_names[cluster_labels == up_cluster]
down_cells = rna.obs_names[cluster_labels == down_cluster]

rxn_de = mern_support.wilcoxon_test(
    enzyme_acts_mean.T,
    up_cells,
    down_cells,
)

rxn_de_info = rxn_info.join(rxn_de)

edge_color_limits = (-2, 2)

edge_values = (
    rxn_de_info["cohens_d"]
    .replace([np.inf, -np.inf], np.nan)
    .fillna(0)
    .clip(lower=edge_color_limits[0], upper=edge_color_limits[1])
)

fig = mern_support.kegg_pathway_plot(
    pathway=pathway_id_to_plot,
    dataset=dataset,
    rna=rna,
    user_labels=edge_values,
    edge_width=4,
    edge_color_limits=edge_color_limits,
    color_scale="RdBu_r",
    show_node_labels=True,
)

fig.update_layout(
    title=(
        f"Glycolysis / Gluconeogenesis: Reaction Activity "
        f"Cluster {up_cluster} vs Cluster {down_cluster}"
    )
)

fig.write_image(f"{base_path_figures}/glycolysis_cohensd_plot.png", width=1200, height=800, scale=2)
fig.show()

cluster 4 (highest mean glycolysis activity) vs cluster 1 (lowest), of 5 clusters with at least 100 cells

Glycolysis Cohen's D plot

Identifying DDPs Connected to Glycolysis#

A KEGG pathway shows reactions within a canonical pathway boundary. DDPs can extend this view by linking pathway reactions to other reactions that are coordinated with them across cells.

Here we identify all DDPs that contain at least one reaction annotated to glycolysis/gluconeogenesis. In the next plot, we expand beyond the KEGG glycolysis pathway by including all reactions from those glycolysis-associated DDPs.

pathway_id_to_plot = "rn00010"

glycolysis_rxns = set(pathway_rxns[pathway_id_to_plot])

glycolysis_ddps = sorted(
    rxn_to_ddp.loc[
        rxn_to_ddp.index.intersection(glycolysis_rxns)
    ]
    .dropna()
    .unique(),
    key=lambda x: int(str(x).split("_")[-1]) if str(x).split("_")[-1].isdigit() else str(x),
)

print("Glycolysis-associated DDPs:", glycolysis_ddps)

Glycolysis-associated DDPs: ['ddp_2', 'ddp_8', 'ddp_56', 'ddp_61', 'ddp_81', 'ddp_86', 'ddp_92', 'ddp_96', 'ddp_97', 'ddp_108']

expanded_rxns = sorted(
    glycolysis_rxns | set().union(*[set(ddp_rxns[ddp]) for ddp in glycolysis_ddps])
)

print("Glycolysis reactions:", len(glycolysis_rxns))
print("Expanded DDP reactions:", len(expanded_rxns))

Glycolysis reactions: 23
Expanded DDP reactions: 97

ddp_edge_labels = rxn_to_ddp.reindex(expanded_rxns)

with ddp_pathway_palette(ddp_edge_labels):
    fig = mern_support.custom_pathway_plot(
        rxns=expanded_rxns,
        dataset=dataset,
        rna=rna,
        user_labels=ddp_edge_labels,
        edge_width=4,
        ddp_labels=rxn_to_ddp,
        ddp_color_map=ddp_color_map_for(ddp_edge_labels),
        ddp_style="soft_contours",
        show_node_labels=True,
        node_label_fraction=0.2,
    )

fig.update_layout(
    title="Expanded Glycolysis-Associated DDP Network"
)

fig.write_image(f"{base_path_figures}/glycolysis_ddp_expanded_plot.png", width=1200, height=800, scale=2)
fig.show()

Glycolysis DDP Expanded plot

Comparing Metabolic Coordination in a Folr1 Knockout#

Everything above describes wild-type (WT) embryogenesis. We now bring in a second dataset: embryos carrying a Folr1 knockout, which severely restricts the embryo’s ability to capture dietary folate.

The metabolic clustering computed above resolves five metabolic states, one of which is almost entirely depleted in the knockout. A whole-dataset comparison would therefore be dominated by that shift in composition, so we compare a single metabolic state across genotypes instead.

Within a matched state, the interesting question is not only which reactions change in mean activity, but which reactions stop moving together. A DDP is defined by coordination, so a perturbation can leave average activity almost untouched while breaking the correlation structure that held a DDP together. The weakest-link analysis is the tool for that comparison, applied here across genotypes rather than across metabolic clusters.

The comparison needs two components:

  1. A shared cell population. The knockout cells are mapped onto the wild-type metabolic states computed above, so we can compare the same metabolic state across genotypes instead of comparing whole datasets.

  2. Replicate-averaged reaction activity for each genotype, decoded from its own MeRN replicates with the analysis behind the paper’s Figure 3.

Loading Reaction Activity for Both Genotypes#

Each genotype’s reaction activity is the replicate ensemble of the models trained on that genotype.

  • WT: enzyme_acts_mean, decoded above from the 5 replicates included with this tutorial (8 decodes each), to demonstrate averaging reaction activity over repeated decodes and over separately trained replicates.

  • Folr1: the knockout models are not included, so its activity is loaded from a precomputed table. Following the paper’s Methods, it is the ensemble mean over 8 Monte Carlo decodes for each of the 10 knockout replicates, restricted to the subsampled knockout cells.

# WT: the replicate-averaged activity decoded above with average_enzyme_activity.
wt_cells = rna.obs_names
wt_enzyme_acts_mean = enzyme_acts_mean.loc[wt_cells]

# Folr1: the knockout's own replicate ensemble, decoded from the models trained on the knockout.
folr1_enzyme_acts_mean = pd.read_csv(
    f'{base_path_data}/{folr1_activity_file}', index_col=0
)
folr1_cells = folr1_enzyme_acts_mean.index

assert not folr1_cells.isin(wt_cells).any(), 'the two genotype tables overlap'

print('WT activity:   ', wt_enzyme_acts_mean.shape)
print('Folr1 activity:', folr1_enzyme_acts_mean.shape)

WT activity:    (6000, 1213)
Folr1 activity: (5000, 1213)

Mapping the Knockout Cells onto the Same Metabolic States#

The wild-type clustering computed above is the reference and is never recomputed.

This follows the procedure of Dias et al.: every Folr1 cell is embedded with the wild-type model, placing both genotypes in one metabolic latent space; each knockout cell then takes the majority label of its 10 closest wild-type neighbours among its 500 nearest neighbours in that shared space. A cell is labelled not_majority when no state holds at least 5 of the 10 votes, and not_enough_neighbors when fewer than 10 of those 500 neighbours are wild-type.

from sklearn.neighbors import NearestNeighbors

n_pool = 500
n_wt_vote = 10
min_votes = 5

# Reuse the WT tutorial-state labels used to construct the DDPs above.

folr1_rna = ad.read_h5ad(f'{base_path_data}/{folr1_adata_file}')[:, features_to_use].copy()
assert list(folr1_rna.var_names) == list(features_to_use)
assert folr1_rna.obs_names.equals(folr1_cells), 'activity table and object disagree on cells'

folr1_latent, _, _ = mern.get_latent_representation(folr1_rna)
print(f'Folr1 metabolic latent: {folr1_latent.shape[0]:,} cells x {folr1_latent.shape[1]} dims')

joint_latent = np.vstack([rna.obsm[metabolic_key], folr1_latent])
is_wt = np.r_[np.ones(rna.n_obs, bool), np.zeros(folr1_rna.n_obs, bool)]

neighbors = NearestNeighbors(n_neighbors=n_pool).fit(joint_latent)
_, neighbor_idx = neighbors.kneighbors(folr1_latent)

# For each knockout cell, take the 10 closest wild-type cells among its 500 neighbours and assign
# the most common state if it holds at least 5 of the 10 votes
wt_label_array = wt_cluster_labels.to_numpy()
assignments = []
for row in neighbor_idx:
    wt_neighbors = row[is_wt[row]][:n_wt_vote]
    if len(wt_neighbors) < n_wt_vote:
        assignments.append('not_enough_neighbors')
        continue
    # np.unique returns labels sorted, so argmax breaks ties toward the first label
    states, counts = np.unique(wt_label_array[wt_neighbors], return_counts=True)
    winner = counts.argmax()
    assignments.append(states[winner] if counts[winner] >= min_votes else 'not_majority')

folr1_cluster_labels = pd.Series(assignments, index=folr1_rna.obs_names)

# A UMAP of the shared embedding
joint = ad.AnnData(obs=pd.DataFrame(index=list(rna.obs_names) + list(folr1_rna.obs_names)))
joint.obsm[metabolic_key] = joint_latent
sc.pp.neighbors(joint, use_rep=metabolic_key)
sc.tl.umap(joint)
wt_umap = joint.obsm['X_umap'][:rna.n_obs]
folr1_umap = joint.obsm['X_umap'][rna.n_obs:]

state_counts = pd.DataFrame({
    'WT': wt_cluster_labels.value_counts(),
    'Folr1': folr1_cluster_labels.value_counts(),
}).fillna(0).astype(int).sort_index()
state_counts.index.name = 'metabolic state'
display(state_counts)

INFO     Input AnnData not setup with scvi-tools. attempting to transfer AnnData setup                             
Folr1 metabolic latent: 5,000 cells x 25 dims
WT Folr1
metabolic state
0 2545 1239
1 1072 1545
2 1040 859
3 754 1102
4 589 25
not_majority 0 230 
reset_plotting_defaults()

states = sorted(set(wt_cluster_labels) | set(folr1_cluster_labels))
state_colors = dict(zip(states, sns.color_palette('tab10', n_colors=len(states)).as_hex()))

wt_germ = rna.obs['germlayer'].astype(str)
folr1_germ = folr1_rna.obs['germlayer'].astype(str)
germlayers = [g for g in GERMLAYER_PALETTE if g in set(wt_germ) | set(folr1_germ)]
missing = (set(wt_germ) | set(folr1_germ)) - set(germlayers)
assert not missing, f'germ layers without a palette color: {missing}'

rows = [
    ('metabolic state', states, state_colors, wt_cluster_labels, folr1_cluster_labels),
    ('germ layer', germlayers, GERMLAYER_PALETTE, wt_germ, folr1_germ),
]

fig, axes = plt.subplots(
    2, 2, figsize=(11, 9.5), sharex=True, sharey=True, layout='constrained'
)
for row_i, (row_axes, (legend_title, categories, palette, wt_labels, folr1_labels)) in enumerate(
    zip(axes, rows)
):
    panels = [
        ('WT', wt_umap, wt_labels),
        ('Folr1 KO', folr1_umap, folr1_labels),
    ]
    for ax, (title, coords, labels_here) in zip(row_axes, panels):
        for category in categories:
            in_category = (labels_here == category).to_numpy()
            if not in_category.any():
                continue
            ax.scatter(
                coords[in_category, 0],
                coords[in_category, 1],
                s=1.5,
                alpha=0.5,
                color=palette[category],
                label=category if ax is row_axes[0] else None,
                rasterized=True,
            )
        if row_i == 0:
            ax.set_title(title)
        if row_i == len(rows) - 1:
            ax.set_xlabel('metabolic UMAP 1')
        ax.set_xticks([])
        ax.set_yticks([])
    row_axes[0].set_ylabel('metabolic UMAP 2')
    handles, labels = row_axes[0].get_legend_handles_labels()
    row_axes[1].legend(
        handles, labels, title=legend_title, markerscale=8,
        bbox_to_anchor=(1.02, 0.5), loc='center left', ncol=1,
    )
sns.despine(fig, left=True, bottom=True)
plt.show()
../_images/b0509f15f03c60ec32d0fc9f754a5f0226c4d71e2ec1cbdfc31d5c9e4fd1d9d3.png

Reusing WT DDPs for the Matched-State Comparison#

We reuse the epiblast-like tutorial state selected above, together with its WT DDPs, correlation matrix, and linkage. Only the knockout correlation matrix and hierarchy are calculated here, using cells mapped onto that same WT state, the same reaction order, and the same DDP threshold.

# Reuse the WT result; do not recluster it for the weakest-link analysis.
wt_ddps = rxn_to_ddp
wt_ddp_rxns = ddp_rxns

folr1_state_cells = folr1_cluster_labels.index[folr1_cluster_labels == comparison_state]
print(f'comparison state {comparison_state}: {len(wt_state_cells)} WT cells, '
      f'{len(folr1_state_cells)} Folr1 cells')

folr1_corr = folr1_enzyme_acts_mean.loc[
    folr1_state_cells, wt_corr.columns
].corr(method='spearman')
folr1_ddps, folr1_ddp_rxns, _, folr1_linkage = mern_support.calculate_ddps(
    met_g, folr1_corr, min_corr=ddp_min_corr, min_size=ddp_min_size
)

print(f'WT DDPs in state {comparison_state}:    {len(wt_ddp_rxns)}')
print(f'Folr1 DDPs in state {comparison_state}: {len(folr1_ddp_rxns)}')

folate_rxns = ['rn:R04325 rn:R06974', 'rn:R04209']
folate_ddps = sorted(set(wt_ddps.reindex(folate_rxns).dropna()))
print('folate reactions:', {r: wt_ddps.get(r) for r in folate_rxns})
folate_ddp = folate_ddps[0] if len(folate_ddps) == 1 else None
if folate_ddp is None:
    print('the two folate reactions are not in a single WT DDP:', folate_ddps)

comparison state 2: 1040 WT cells, 859 Folr1 cells
WT DDPs in state 2:    115
Folr1 DDPs in state 2: 140
folate reactions: {'rn:R04325 rn:R06974': 'ddp_89', 'rn:R04209': 'ddp_89'}

WT DDPs Across Central Carbon and Nucleotide Metabolism#

Before scoring what breaks, it helps to see what the WT programs are. We take every DDP that touches glycolysis, pentose phosphate, purine, or pyrimidine metabolism and draw their compound connections, colored by DDP membership. Each DDP is a connected, coordinated block. The plotting function labels a subset of metabolites across the full layout to limit overlap; hover over any node to see its name.

folate_ddp sits in purine metabolism: its reactions are consecutive steps of de novo purine synthesis, the route that consumes a pentose phosphate pathway product to build the purine ring.

Only WT DDPs appear here. The knockout’s own DDPs are never used as groups: the weakest-link analysis below keeps the WT DDPs fixed and uses the knockout’s correlations and clustering hierarchy only to score how well each WT DDP holds together.

pathways_of_interest = ['rn00010', 'rn00030', 'rn00230', 'rn00240']

keep_ddps = set()
keep_rxns = set()
for rxn in wt_corr.index:
    if rxn not in rxn_info.index:
        continue
    rxn_pathways = str(rxn_info.loc[rxn, 'Pathways']).split(';')
    if not set(rxn_pathways) & set(pathways_of_interest):
        continue
    keep_rxns.add(rxn)
    if rxn in wt_ddps.index and wt_ddps[rxn] is not None and not pd.isnull(wt_ddps[rxn]):
        keep_ddps.add(wt_ddps[rxn])

for ddp in keep_ddps:
    keep_rxns.update(wt_ddp_rxns[ddp])

keep_rxns = sorted(keep_rxns)
print(f'{len(keep_ddps)} DDPs touch the selected pathways, spanning {len(keep_rxns)} reactions')

24 DDPs touch the selected pathways, spanning 258 reactions

wt_labels = wt_ddps.where(wt_ddps.index.isin(keep_rxns), None)

with ddp_pathway_palette(wt_labels):
    fig = mern_support.custom_pathway_plot(
        rxns=keep_rxns,
        dataset=dataset,
        rna=rna,
        user_labels=wt_labels,
        show_node_labels=True,
        node_label_fraction=0.2,
        node_label_font_size=9,
    )
fig.update_layout(title='Selected reactions colored by WT DDP membership')
fig.write_image(f"{base_path_figures}/folr1_ddp_plot.png", width=1200, height=800, scale=2)
fig.show()

Folr1 DDP Plot

Contents