Visualizing Inferred GRN#

RegDiffusion hands you back an adjacency matrix. On a 5,000-gene run that is 25 million numbers, and no amount of squinting at a heatmap will tell you whether the network learned anything biological. The fastest way to build (or lose) confidence in an inferred network is to look at it and ask whether genes that belong together ended up together.

This tutorial walks through that end to end on a real dataset: pull healthy human PBMCs from CELLxGENE, infer a GRN over 5,000 genes, annotate the genes with Reactome pathways, and render the result with lightgraph.

Note

Network plots in regdiffusion are drawn by lightgraph, which renders on an HTML canvas with viewport culling. This matters more than it sounds: the older pyvis/vis.js backend stalls for the better part of a minute before the first frame at 5,000 nodes, which effectively ruled out looking at a whole network at once. lightgraph stays interactive at that size, so the whole-network view below is something you can actually pan around in.

Step 1. Get the data#

We use the NYGC multimodal PBMC dataset from Hao et al., Cell 2021 — the Seurat v4 CITE-seq reference, 161,764 healthy human PBMCs, deeply annotated, and a de facto reference for immune cell types. It is available through the CELLxGENE Census API:

>>> import cellxgene_census
>>> import scanpy as sc
>>> import numpy as np
>>>
>>> DATASET_ID = "ed5d841d-6346-47d4-ab2f-7119ad7e3a35"
>>>
>>> with cellxgene_census.open_soma(census_version="2023-12-15") as census:
>>>     adata = cellxgene_census.get_anndata(
>>>         census,
>>>         organism="Homo sapiens",
>>>         obs_value_filter=f"dataset_id == '{DATASET_ID}' and is_primary_data == True",
>>>         obs_column_names=["cell_type", "assay", "donor_id", "disease"],
>>>     )
>>> adata.shape
(161764, 60664)

Census ships raw counts indexed by Ensembl ID, so we switch to gene symbols (which is what pathway databases speak) and trim the gene families that soak up variance without carrying much regulatory signal:

>>> adata.var["ensembl_id"] = adata.var_names
>>> adata.var_names = adata.var["feature_name"].astype(str)
>>> adata.var_names_make_unique()
>>>
>>> mask = ~(
>>>     adata.var_names.str.startswith("MT-")
>>>     | adata.var_names.str.startswith(("RPL", "RPS"))
>>>     | adata.var_names.str.match(r"^HB[ABDEGQZ]\d?$")
>>> )
>>> adata = adata[:, mask].copy()
>>> sc.pp.filter_genes(adata, min_cells=int(0.01 * adata.n_obs))
>>>
>>> sc.pp.normalize_total(adata, target_sum=1e4)
>>> sc.pp.log1p(adata)
>>> sc.pp.highly_variable_genes(adata, n_top_genes=5000, flavor="seurat")
>>> adata = adata[:, adata.var.highly_variable].copy()
>>> adata.shape
(161764, 5000)

Step 2. Infer the network#

Nothing special here — this is the standard RegDiffusion workflow. Passing cell_types lets the model condition on the 30 annotated cell types in this dataset.

>>> import regdiffusion as rd
>>>
>>> X = np.asarray(adata.X.todense(), dtype=np.float32)
>>> codes = adata.obs["cell_type"].astype("category").cat.codes.values.astype(np.int64)
>>>
>>> trainer = rd.RegDiffusionTrainer(X, cell_types=codes, n_steps=1000, seed=123)
>>> trainer.train()
>>> grn = trainer.get_grn(adata.var_names.to_numpy())
>>> grn.to_hdf5("pbmc_nygc_grn.hdf5", as_sparse=True)

On an RTX PRO 6000 the 1,000 training steps take about 13 seconds; loading and densifying the expression matrix takes considerably longer than the training itself.

A first sanity check before any plotting — who are the hub genes?

>>> from lightgraph import degree
>>>
>>> edges = grn.extract_edgelist(k=2)
>>> edges = edges.assign(weight=edges.weight.abs())
>>> deg = degree(edges[["source", "target", "weight"]].values.tolist())
>>> sorted(deg.items(), key=lambda kv: -kv[1])[:12]
[('LTB', 306.0), ('IL7R', 292.0), ('ETS1', 203.0), ('CD3D', 181.0),
 ('NKG7', 149.0), ('CD3E', 145.0), ('LDHB', 144.0), ('CD79A', 138.0),
 ('EEF1B2', 138.0), ('MS4A1', 131.0), ('CST7', 129.0), ('EEF1A1', 115.0)]

That is a reassuring list: T cell identity (CD3D, CD3E, IL7R, LTB), B cell identity (CD79A, MS4A1), cytotoxic effectors (NKG7, CST7), and the lymphoid transcription factor ETS1. The network put PBMC lineage markers at the center without being told anything about lineages.

Step 3. Group genes by function#

Topology is much easier to judge when nodes carry an independent functional label. We use Reactome’s pathway hierarchy: every gene is annotated to one or more pathways, and each pathway rolls up to one of ~30 top-level pathways like Immune System or Cell Cycle. Rolling annotations up to that top level gives a small, readable palette.

Two files from Reactome are all we need — the gene sets and the parent/child relations:

wget https://reactome.org/download/current/ReactomePathways.gmt.zip
wget https://reactome.org/download/current/ReactomePathwaysRelation.txt
unzip ReactomePathways.gmt.zip

>>> from collections import defaultdict, Counter
>>>
>>> name, gene2path = {}, defaultdict(set)
>>> for line in open("ReactomePathways.gmt"):
>>>     f = line.rstrip("\n").split("\t")
>>>     name[f[1]] = f[0]
>>>     for g in f[2:]:
>>>         gene2path[g].add(f[1])
>>>
>>> parent, children = defaultdict(set), defaultdict(set)
>>> for line in open("ReactomePathwaysRelation.txt"):
>>>     p, c = line.rstrip("\n").split("\t")
>>>     if p.startswith("R-HSA") and c.startswith("R-HSA"):
>>>         parent[c].add(p)
>>>         children[p].add(c)
>>>
>>> roots = {p for p in set(parent) | set(children) | set(name) if not parent[p]}
>>>
>>> def ancestors(pid, cache={}):
>>>     if pid not in cache:
>>>         cache[pid] = ({pid} if pid in roots or not parent[pid]
>>>                       else set().union(*(ancestors(p) for p in parent[pid])))
>>>     return cache[pid]

A gene often reaches several top-level pathways. We let the roots vote — the one covering most of the gene’s pathways wins — and break ties toward the smaller root, so genes do not all avalanche into Metabolism and Signal Transduction:

>>> root_size = Counter()
>>> for g, pids in gene2path.items():
>>>     for r in set().union(*(ancestors(p) for p in pids)):
>>>         root_size[r] += 1
>>>
>>> def assign(symbols):
>>>     out = {}
>>>     for s in symbols:
>>>         votes = Counter()
>>>         for p in gene2path.get(s, ()):
>>>             for a in ancestors(p):
>>>                 votes[a] += 1
>>>         if votes:
>>>             best = min(votes, key=lambda r: (-votes[r], root_size[r], r))
>>>             out[s] = name.get(best, best)
>>>     return out
>>>
>>> groups_all = assign(grn.gene_names)
>>> len(groups_all), len(grn.gene_names)
(2678, 5000)

Roughly 54% of our genes get a Reactome annotation. That is normal and worth being upfront about: highly variable genes in single cell data are full of lncRNAs, immunoglobulin and TCR segments, and poorly characterized ORFs that no curated database covers. Keep the twelve largest groups and simply leave everything else out of the dictionary — genes missing from node_group_dict render in a neutral gray and stay out of the legend, which is much more honest than dumping half the network into a colored bucket called “Other”:

>>> top12 = [g for g, _ in Counter(groups_all.values()).most_common(12)]
>>> groups = {g: p for g, p in groups_all.items() if p in top12}
>>> Counter(groups.values()).most_common()
[('Immune System', 403), ('Signal Transduction', 366),
 ('Gene expression (Transcription)', 338), ('Metabolism', 335),
 ('Metabolism of proteins', 179), ('Disease', 154), ('Cell Cycle', 142),
 ('Metabolism of RNA', 125), ('Vesicle-mediated transport', 76),
 ('Transport of small molecules', 70), ('DNA Repair', 66),
 ('Developmental Biology', 64)]

Give the pathways fixed, meaningful colors#

This tutorial shows several views of the same network — the whole graph, a local neighborhood, and the same analysis rerun on a different cell type. By default lightgraph assigns palette colors in sorted-group order, so a group that is absent from one view shifts the colors of everything after it and the figures stop being comparable.

Pin the colors instead. Assigning hues by biological family also means related pathways read as related — the two genome-maintenance groups are both violet, the three metabolism groups are all green, and so on:

>>> PATHWAY_COLORS = {
>>>     # defense and pathology
>>>     'Immune System':                   '#e34948',
>>>     'Disease':                         '#ff756f',
>>>     # stress / alert
>>>     'Cellular responses to stimuli':   '#eda100',
>>>     # genome maintenance
>>>     'Cell Cycle':                      '#4a3aa7',
>>>     'DNA Repair':                      '#918bff',
>>>     # metabolism
>>>     'Metabolism':                      '#008300',
>>>     'Metabolism of proteins':          '#36a231',
>>>     'Metabolism of RNA':               '#56bf50',
>>>     # information flow
>>>     'Gene expression (Transcription)': '#2a78d6',
>>>     # signaling
>>>     'Signal Transduction':             '#e87ba4',
>>>     # trafficking
>>>     'Vesicle-mediated transport':      '#eb6834',
>>>     'Transport of small molecules':    '#ff8552',
>>>     # development
>>>     'Developmental Biology':           '#1baf7a',
>>> }
>>> GROUP_ORDER = list(PATHWAY_COLORS)   # union across every figure below

Pass both to every figure. group_order reserves a palette slot for each group even when that group is missing from the current subset, so nothing shifts:

>>> COMMON = dict(group_colors=PATHWAY_COLORS, group_order=GROUP_ORDER)

Note

Thirteen categories is more than color alone can safely carry. Under red-green color vision deficiency some of these pairs are genuinely indistinguishable, and no palette fixes that — it is a limit of the channel, not of the color choices. The design leans on that deliberately: confusions are pushed inside a family, where they cost least (mistaking Metabolism of RNA for Metabolism of proteins is far less damaging than mistaking either for Cell Cycle). Identity is never carried by color alone in these figures — every group is named with its size in the legend, clicking a legend entry isolates that group, hovering a node names its group, and the search box finds genes directly.

Step 4. Look at the whole network#

Now we can put all 5,000 genes on screen at once. We keep the top 2 edges per gene (10,000 edges), color by Reactome group, and size each node by its centrality so the regulatory hubs are visible at a glance.

PageRank is the natural centrality here — cheap to compute even at this size, and it agrees with the degree ranking on the top hubs. One catch: PageRank on this graph spans roughly 375× from the smallest gene to the largest, with a very long tail. net_vis min-max normalizes the metric linearly, so feeding the raw values would leave the median gene at 0.002 of the size range — i.e. everything minimum-size except a handful of hubs. Taking the log first spreads the middle of the distribution across the size channel:

>>> import numpy as np
>>> from lightgraph import net_vis, pagerank
>>>
>>> pr = pagerank(edges[["source", "target", "weight"]].values.tolist())
>>> centrality = {g: float(np.log10(v)) for g, v in pr.items()}
>>>
>>> v = net_vis(
>>>     edges=edges,
>>>     node_groups=groups,
>>>     node_metric=centrality,
>>>     metric_map="size",
>>>     metric_label="PageRank (log)",
>>>     metric_size_range=(5, 32),
>>>     node_color="#e0e0e0",            # unannotated genes
>>>     edge_weight_to_opacity=True,
>>>     weight_opacity_range=(0.008, 0.10),
>>>     simulation_strength=12000,
>>>     link_distance=45,
>>>     warmup_ticks=600,
>>>     show_labels=False,
>>>     show_ellipses=False,
>>>     show_statistics=True,
>>>     height="720px",
>>>     **COMMON,                        # pinned pathway colors
>>> )
>>> v.save("grn_pbmc_full.html")

With 10,000 edges behind 5,000 sizeable nodes, the edge layer has to get out of the way — hence the very low weight_opacity_range. The edges still do their job of holding the layout together and darkening where the graph is dense, without turning the picture into a gray wash.

The top genes by PageRank are the same lineage markers the degree ranking found: IL7R, LTB, ETS1, CD3D, NKG7, CD3E, CST7, CD79A, MS4A1.

Warning

Note what we are not doing: lightgraph has a group_attraction parameter that pulls same-group nodes toward a shared centroid. Turning it up would make the pathway groups look beautifully separated — but it would be manufacturing that separation from the labels rather than revealing it from the edges. We leave it at its default so that any clustering you see comes from the inferred network alone.

Things worth trying in the widget above: hover a node to fade everything outside its immediate neighborhood, double-click a node to isolate its neighborhood entirely (double-click empty space or hit Escape to come back), and click a legend entry to focus one pathway group. Labels are off by default at this size — zoom in and they appear.

Does the picture actually mean anything?#

Eyeballing a hairball is a good way to fool yourself, so it is worth checking the impression numerically. If functionally related genes really are wired together, edges should connect same-group genes more often than chance. Shuffle the pathway labels among annotated genes and compare:

>>> both = edges[edges.source.isin(groups) & edges.target.isin(groups)]
>>> sg, tg = both.source.map(groups), both.target.map(groups)
>>> observed = (sg == tg).mean()
>>>
>>> rng = np.random.default_rng(0)
>>> keys, vals = list(groups), list(groups.values())
>>> null = np.mean([
>>>     (both.source.map(m) == both.target.map(m)).mean()
>>>     for m in (dict(zip(keys, rng.permutation(vals))) for _ in range(200))
>>> ])
>>> observed, null, observed / null
(0.197, 0.118, 1.67)

Same-group edges are 1.67× more common than chance (z ≈ 10). Broken down by group — observed vs. its own shuffled null — the effect is very unevenly distributed:

Reactome group	within-group	null	ratio
DNA Repair	0.052	0.013	3.95
Cell Cycle	0.094	0.030	3.10
Immune System	0.217	0.094	2.32
Vesicle-mediated transport	0.024	0.014	1.63
Disease	0.053	0.034	1.59
Signal Transduction	0.127	0.085	1.49
Transport of small molecules	0.019	0.014	1.38
Metabolism of RNA	0.023	0.026	0.89
Metabolism of proteins	0.035	0.040	0.87
Metabolism	0.066	0.077	0.86
Gene expression (Transcription)	0.064	0.079	0.80
Developmental Biology	0.007	0.012	0.55

The coherent groups are exactly the ones that should be coherent in resting PBMCs: DNA repair and cell cycle are tightly co-regulated programs, and immune function is the biology this tissue is built around. The groups at or below 1.0 are the broad grab-bag categories — Metabolism and Gene expression (Transcription) each cover hundreds of genes spanning unrelated processes, so there is no reason for them to occupy one region of the graph, and they do not.

In other words, the network is not uniformly “functionally organized”. It is organized around specific, tightly co-regulated programs, and the visualization is showing you something real rather than a layout artifact.

Should you filter out weak edges?#

Worth being clear about what extract_edgelist(k=2) does: it is a per-gene rank filter, not a strength filter. Every gene contributes its 2 best edges whether those edges are strong or barely there. Since the weight distribution is heavily skewed, that lets a lot of weak links in:

>>> np.quantile(edges.weight, [0, 0.25, 0.5, 0.75, 0.9, 1])
array([ 1.136,  1.545,  2.102,  3.988,  7.578, 15.32 ])

Adding a weight threshold on top sharpens the functional signal considerably — here is the same-group enrichment from above, recomputed at increasing cutoffs:

cutoff	edges	genes kept	enrichment
none	10,000	5,000	1.67×
25th pct (1.55)	7,504	3,921	1.78×
median (2.10)	5,007	2,602	2.00×
75th pct (3.99)	2,502	1,291	2.64×
90th pct (7.58)	1,001	525	3.23×

It is tempting to read that column as “strong edges are simply more trustworthy” and threshold aggressively. Don’t — look at which genes survive before you believe it. Here is how each pathway group is thinned at each cutoff (percentage of the group’s genes still on the plot):

Reactome group	none	median	75th pct	90th pct
Immune System	403	331 (82%)	255 (63%)	145 (36%)
Signal Transduction	366	268 (73%)	164 (45%)	64 (17%)
Disease	154	103 (67%)	56 (36%)	27 (18%)
Developmental Biology	64	44 (69%)	28 (44%)	13 (20%)
Metabolism	335	165 (49%)	74 (22%)	29 (9%)
Gene expression (Transcription)	338	124 (37%)	49 (14%)	17 (5%)
Cell Cycle	142	57 (40%)	22 (15%)	2 (1%)
DNA Repair	66	26 (39%)	4 (6%)	1 (2%)

Weight filtering is not biologically neutral. It preferentially keeps the dominant program of the tissue — in resting PBMCs, immune identity genes carry by far the largest weights — and it deletes Cell Cycle and DNA Repair almost entirely by the 90th percentile. Those are precisely the two groups with the highest functional coherence in the table above (3.10× and 3.95×).

Which means the rising global enrichment is partly an artifact: the graph is collapsing onto one large, internally-consistent immune module, and “same-group edges” rises because one group is swallowing the plot. You are not making the picture more trustworthy so much as making it more about T and B cells.

Reduce `k` instead of thresholding#

There is a better lever, and it is the one extract_edgelist already gives you. Note what k actually does: for each target gene it keeps that gene’s k strongest regulators. It is a rank filter applied per gene, so it is invariant to the scale of any particular gene’s weights — a cell cycle gene keeps its best regulators even if those weights are small in absolute terms. (This is the same primitive as select_top_n_per_target in flashscenic, which builds SCENIC modules on top of RegDiffusion and unions several rank-based selections precisely so that no single global cutoff decides what survives.)

Compare the two ways of getting down to about 5,000 edges:

	edges	genes kept	Immune System	Cell Cycle	DNA Repair	enrichment
`extract_edgelist(k=1)`	5,000	5,000	403	142	66	1.71×
`k=2` + median weight cut	5,007	2,602	331	57	26	2.05×

Same edge budget. The rank-based version keeps every gene and every pathway group completely intact; the weight cut throws away half the genes and 60% of the cell cycle. The weight cut does score a somewhat higher enrichment (2.05× vs 1.71×) — but as established above, a good part of that is the compositional artifact of the immune program crowding out everything else, not a cleaner graph.

So, practical guidance:

To make the plot sparser — lower k. It is scale-free, it costs you nothing in gene coverage, and no program silently vanishes.
To de-emphasize weak links without deleting them — keep every edge and let opacity carry the weight (edge_weight_to_opacity=True, as used above).
To threshold on weight anyway — check your program of interest survives first, and say so in the caption, because the composition of the figure changed and not only its density.

GRN has this built in, applied to the adjacency matrix before you extract anything:

>>> grn.remove_weak_edges(threshold=2.1)
>>> edges_strong = grn.extract_edgelist(k=2)

Note that remove_weak_edges mutates the GRN in place and only ever raises the cutoff, so reload from HDF5 if you want to go back to a weaker threshold.

Step 5. Zoom in on a neighborhood#

Whole-network views are good for gestalt, bad for specifics. For a claim you can check gene by gene, use visualize_local_neighborhood, which pulls the 2.5-hop neighborhood around genes you name. Let us take the two B cell hubs from the degree ranking:

>>> g = grn.visualize_local_neighborhood(
>>>     ["MS4A1", "CD79A"], k=20, node_group_dict=groups, height="640px",
>>>     **COMMON,
>>> )
>>> g.save("grn_pbmc_bcell.html")

The two genes you asked about are drawn at double size, and edge width follows the absolute regulatory strength of each link. The neighborhood comes back with 16 genes:

BANK1, CD22, CD79A, CD79B, HLA-DQA1, IGHD, IGHM, IGKC, LINC00926,
MEF2C, MS4A1, PAX5, RALGPS2, TCL1A, TNFRSF13C, VPREB3

Every one of these is a B cell gene. CD79A/CD79B/CD22/MS4A1 are B cell receptor complex and surface markers; IGHD/IGHM/IGKC are immunoglobulin chains; PAX5 and MEF2C are the master transcription factors of B cell commitment; TCL1A and VPREB3 mark naive and pre-B stages; TNFRSF13C (BAFF-R) drives B cell survival. Not a single T cell, myeloid, or NK gene got pulled in, even though the run had 5,000 genes to choose from and no cell type labels were used to build the neighborhood.

Note also that LINC00926 and RALGPS2 sit right in this module. Neither carries a Reactome annotation, so both render in the neutral gray — but their company is a strong hint, and the literature bears it out. LINC00926 is a lncRNA whose expression is restricted to the B cell lineage (highest in naive B cells, down-regulated on activation), and the Human Protein Atlas classifies RALGPS2 as B cell group enriched among immune cells. This is the practical payoff of an inferred GRN: it places genes your annotation database has nothing to say about right next to the ones it does.

Step 6. A different cell type, a different network#

Everything above was inferred across all 30 PBMC cell types at once, so the strongest signal in the data is lineage identity — which is why the hubs are T and B cell markers. A good way to check that the method is reading biology rather than reproducing one fixed answer is to run it inside a single cell type, where lineage identity is constant and cannot dominate.

CD14+ monocytes are the largest population here (42,690 cells) and about as far from the T/B axis as this dataset goes. The pipeline is identical, with two changes: highly variable genes are re-selected within monocytes, so they capture variation across monocyte states rather than between lineages, and there is no cell type covariate left to condition on.

>>> with cellxgene_census.open_soma(census_version="2023-12-15") as census:
>>>     mono = cellxgene_census.get_anndata(
>>>         census, organism="Homo sapiens",
>>>         obs_value_filter=(f"dataset_id == '{DATASET_ID}' and is_primary_data == True "
>>>                           "and cell_type == 'CD14-positive monocyte'"),
>>>     )
>>> # ...same gene filtering, normalization and HVG selection as Step 1...
>>> trainer = rd.RegDiffusionTrainer(X_mono, n_steps=1000, seed=123)   # no cell_types
>>> trainer.train()
>>> mono_grn = trainer.get_grn(mono.var_names.to_numpy())

The monocyte figure below is rendered with the same **COMMON colors, which is what makes it comparable to the whole-PBMC one: Immune System is the same red in both, Cell Cycle the same violet. Two groups differ between the runs — DNA Repair drops out of the monocyte top twelve and Cellular responses to stimuli enters — and because group_order reserves a slot per group, that swap changes nothing else on screen.

Only 2,739 of the 5,000 genes overlap with the all-PBMC run, so this is genuinely a different gene space. The hubs change completely:

	top genes by PageRank
All PBMC	`IL7R`, `LTB`, `ETS1`, `CD3D`, `NKG7`, `CD3E`, `CST7`, `CD79A`, `MS4A1`
CD14+ monocytes	`S100A8`, `TXNIP`, `S100A9`, `IL1B`, `HLA-DPA1`, `HLA-DPB1`, `MNDA`, `CD74`, `HLA-DRA`, `S100A12`, `VCAN`

Not a single gene in common, and the monocyte list is textbook myeloid biology: the S100A8/S100A9/S100A12 alarmins, IL1B for inflammatory activation, and the MHC class II machinery (HLA-DRA, HLA-DPA1, HLA-DPB1, CD74) that monocytes use to present antigen.

The pathway coherence analysis from Step 4 tells the sharper story. Overall enrichment is essentially unchanged — 1.63× in monocytes vs 1.67× across all PBMCs — so both networks are organized to the same degree. But which programs are coherent is almost entirely different:

Reactome group	monocytes	all PBMC
Cellular responses to stimuli	4.64×	not in top 12
Transport of small molecules	2.71×	1.28×
Immune System	2.16×	2.27×
Metabolism of RNA	1.46×	0.92×
Signal Transduction	1.33×	1.52×
Cell Cycle	1.00×	3.12×
DNA Repair	too few edges	4.12×

Cell cycle coherence collapses from 3.12× to exactly chance, and DNA repair does not even have enough edges left to test. That is the correct answer: circulating monocytes are terminally differentiated and essentially non-proliferating, so there is no coordinated cell cycle program for the model to find. What replaces it at the top is Cellular responses to stimuli at 4.64× — stress and stimulus response, which is precisely what distinguishes one circulating monocyte from another.

The method is not applying a fixed template. Run it on a mixed population with proliferating cells and cell cycle comes out coherent; run it inside a post-mitotic population and that structure correctly disappears.