Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Gene regulatory networks

🧠 Key takeaways

Multiome data allow regions to be linked to genes and TFs to regions in the same cells, but the resulting networks are hypotheses that differ between methods.

Region-gene links combine genomic proximity with the correlation of accessibility and expression, and TF-region links require a binding motif.

Only a small fraction of inferred links is supported by prior knowledge, so networks need validation, ideally by perturbing the TFs.

⚙️ Environment setup
Steps
yml

  1. Install conda:

    • Before creating the environment, ensure that conda is installed on your system.

  2. Save the yml content:

    • Copy the content from the yml tab into a file named environment.yml.

  3. Create the environment:

    • Open a terminal or command prompt.

    • Run the following command:

      conda env create -f environment.yml

  4. Activate the environment:

    • After the environment is created, activate it using:

      conda activate <environment_name>

    • Replace <environment_name> with the name specified in the environment.yml file. In the yml file it will look like this:

      name: <environment_name>

  5. Verify the installation:

    • Check that the environment was created successfully by running:

      conda env list
🗄️ Get data and notebooks

This book uses lamindb to store, share, and load datasets and notebooks using the theislab/sc-best-practices instance. We acknowledge free hosting from Lamin Labs.

  1. Install lamindb

    • Install the lamindb Python package:

    pip install lamindb

  2. Optionally create a lamin account

  3. Verify your setup

    • Run the lamin connect command:

    import lamindb as ln

ln.Artifact.connect("theislab/sc-best-practices").df()

    You should now see up to 100 of the stored datasets.

  4. Accessing datasets (Artifacts)

    • Search for the datasets on the Artifacts page

    • Load an Artifact and the corresponding object:

    import lamindb as ln

af = ln.Artifact.connect("theislab/sc-best-practices").get(
    key="key_of_dataset", is_latest=True
)
obj = af.load()

    The object is now accessible in memory and is ready for analysis. Adapt the lamindb.Artifact.connect("theislab/sc-best-practices").get("SOMEIDXXXX") suffix to get respective versions.

  5. Accessing notebooks (Transforms)

    lamin load <notebook url>

    which will download the notebook to the current working directory. Analogously to Artifacts, you can adapt the suffix ID to get older versions.

Motivation

Transcription factors (TFs) regulate genes by binding to accessible regulatory regions, such as promoters and enhancers, whose accessibility changes with their activity Spitz & Furlong, 2012. Multiome data measure accessibility and expression in the same cells, which allows us to link regions to the genes whose expression follows their accessibility, and TFs to the regions that contain their motifs. Together, these links form a network from TFs through regions to target genes, the idea behind methods such as SCENIC+ Bravo González-Blas et al., 2023.

Such networks are hypotheses rather than measurements: the networks of different methods barely overlap, change with random seeds and preprocessing, and predict the effects of perturbations poorly Mompel et al., 2024. For TF activities from gene expression alone, see the gene regulatory network chapter of the RNA section.

import lamindb as ln
import pandas as pd
import scanpy as sc
import snapatac2 as snap

ln.connect("theislab/sc-best-practices")

ln.track("N64Xsy5YJqH4")
Output
→ connected lamindb: theislab/sc-best-practices

→ loaded Transform('N64Xsy5YJqH40000', key='gene_regulatory_networks_atac.ipynb'), started new Run('Zh98aN728L2pCIMK') at 2026-09-29 15:44:48 UTC

→ notebook imports: lamindb-core==2.10.0 pandas==2.3.3 scanpy==1.12.4 snapatac2==2.10.0

We load the peak matrix from the previous chapter and the gene expression of the same cells, and normalize both with the shifted logarithm. We restrict the network to the cell type-specific peaks and to the 2,000 most variable genes, which keeps it small and focused on regulation that differs between cell types.

peak_matrix = ln.Artifact.get(key="chromatin_accessibility/pbmc10k_peaks.h5ad").load()
rna = ln.Artifact.get(
    key="introduction/advanced_data_structures_and_frameworks_muon.h5mu"
).load()["rna"]
rna = rna[peak_matrix.obs_names].copy()
rna.var_names_make_unique()
for data in (peak_matrix, rna):
    sc.pp.normalize_total(data)
    sc.pp.log1p(data)
sc.pp.highly_variable_genes(rna, n_top_genes=2000)
markers = snap.tl.marker_regions(peak_matrix, groupby="cell_type", pvalue=0.01)
regions = sorted({region for peaks in markers.values() for region in peaks})

Linking regions to genes

We connect every region to the genes whose transcription start site lies within 250 kb and keep the links in which accessibility and expression are correlated across cells.

annotation = ln.Artifact.get(
    key="chromatin_accessibility/gencode_v41_GRCh38.gff3.gz"
).cache()
network = snap.tl.init_network_from_annotation(regions, annotation)
variable_genes = set(rna.var_names[rna.var["highly_variable"]])
network = snap.tl.prune_network(
    network, node_filter=lambda node: node.type != "gene" or node.id in variable_genes
)
snap.tl.add_cor_scores(network, gene_mat=rna, peak_mat=peak_matrix)
network = snap.tl.prune_network(
    network, edge_filter=lambda source, target, link: link.cor_score > 0.1
)
Output
2026-09-29 17:45:03 - INFO - Removed 22374 isolated nodes.

  0%|          | 0/1204 [00:00<?, ?it/s]
  4%|▎         | 44/1204 [00:00<00:02, 432.77it/s]
  7%|▋         | 89/1204 [00:00<00:02, 436.39it/s]
 11%|█         | 133/1204 [00:00<00:02, 412.37it/s]
 15%|█▍        | 178/1204 [00:00<00:02, 425.95it/s]
 18%|█▊        | 222/1204 [00:00<00:02, 428.96it/s]
 22%|██▏       | 265/1204 [00:00<00:02, 428.19it/s]
 26%|██▌       | 310/1204 [00:00<00:02, 431.96it/s]
 29%|██▉       | 354/1204 [00:00<00:02, 415.27it/s]
 33%|███▎      | 400/1204 [00:00<00:01, 428.47it/s]
 37%|███▋      | 444/1204 [00:01<00:01, 426.45it/s]
 41%|████      | 489/1204 [00:01<00:01, 430.75it/s]
 44%|████▍     | 533/1204 [00:01<00:01, 428.78it/s]
 48%|████▊     | 576/1204 [00:01<00:01, 426.61it/s]
 52%|█████▏    | 622/1204 [00:01<00:01, 431.36it/s]
 55%|█████▌    | 667/1204 [00:01<00:01, 434.01it/s]
 59%|█████▉    | 711/1204 [00:01<00:01, 419.47it/s]
 63%|██████▎   | 755/1204 [00:01<00:01, 424.42it/s]
 66%|██████▋   | 798/1204 [00:01<00:00, 424.21it/s]
 70%|██████▉   | 841/1204 [00:01<00:00, 424.86it/s]
 74%|███████▍  | 890/1204 [00:02<00:00, 440.00it/s]
 78%|███████▊  | 935/1204 [00:02<00:00, 422.31it/s]
 81%|████████  | 978/1204 [00:02<00:00, 401.90it/s]
 85%|████████▍ | 1021/1204 [00:02<00:00, 408.67it/s]
 88%|████████▊ | 1064/1204 [00:02<00:00, 413.80it/s]
 92%|█████████▏| 1108/1204 [00:02<00:00, 420.65it/s]
 96%|█████████▌| 1152/1204 [00:02<00:00, 425.19it/s]
 99%|█████████▉| 1195/1204 [00:02<00:00, 426.13it/s]
100%|██████████| 1204/1204 [00:02<00:00, 424.78it/s]

2026-09-29 17:45:06 - INFO - Removed 10238 isolated nodes.

f"{network.num_edges()} region-gene links"
'541 region-gene links'

Adding transcription factors

Next, we connect TFs to the linked regions that contain their CIS-BP motif Weirauch et al., 2014 and derive TF-gene links through these regions Zhang et al., 2024. A motif only shows that a TF could bind, so we additionally require that the TF and its target are co-expressed, and we remove links of TFs to their own gene, whose correlation is trivially 1.

fasta = ln.Artifact.get(key="chromatin_accessibility/gencode_v41_GRCh38.fa").cache()
motif_file = ln.Artifact.get(key="chromatin_accessibility/cisBP_human.meme").cache()
cis_bp = {}
for motif in snap.read_motifs(str(motif_file)):
    motif.name = motif.id.split("+")[0]
    if (
        motif.name not in cis_bp
        or motif.info_content() > cis_bp[motif.name].info_content()
    ):
        cis_bp[motif.name] = motif
snap.tl.add_tf_binding(network, motifs=list(cis_bp.values()), genome_fasta=fasta)
grn = snap.tl.link_tf_to_gene(network)
snap.tl.add_cor_scores(grn, gene_mat=rna, peak_mat=peak_matrix)
grn = snap.tl.prune_network(
    grn,
    edge_filter=lambda source, target, link: source != target and link.cor_score > 0.1,
)
edges = pd.DataFrame(
    [
        (grn[source].id, grn[target].id, grn.get_edge_data(source, target).cor_score)
        for source, target in grn.edge_list()
    ],
    columns=["tf", "target", "correlation"],
)
Output
2026-09-29 17:45:07 - INFO - Fetching 524 sequences ...

2026-09-29 17:45:07 - INFO - Searching for the binding sites of 1165 motifs ...

  0%|          | 0/1165 [00:00<?, ?it/s]
  1%|          | 7/1165 [00:00<00:18, 62.18it/s]
  2%|▏         | 24/1165 [00:00<00:09, 121.47it/s]
  3%|▎         | 37/1165 [00:00<00:12, 92.85it/s] 
  4%|▍         | 48/1165 [00:00<00:12, 92.22it/s]
  5%|▌         | 64/1165 [00:00<00:09, 110.39it/s]
  7%|▋         | 76/1165 [00:00<00:09, 111.01it/s]
  8%|▊         | 88/1165 [00:00<00:10, 101.16it/s]
  9%|▉         | 103/1165 [00:00<00:09, 113.91it/s]
 10%|▉         | 115/1165 [00:01<00:09, 112.72it/s]
 11%|█         | 127/1165 [00:01<00:10, 100.78it/s]
 12%|█▏        | 138/1165 [00:01<00:11, 92.49it/s] 
 13%|█▎        | 151/1165 [00:01<00:10, 97.94it/s]
 14%|█▍        | 162/1165 [00:01<00:10, 93.99it/s]
 15%|█▍        | 172/1165 [00:01<00:11, 89.16it/s]
 16%|█▌        | 182/1165 [00:01<00:11, 83.67it/s]
 16%|█▋        | 191/1165 [00:02<00:13, 72.41it/s]
 17%|█▋        | 199/1165 [00:02<00:14, 67.03it/s]
 18%|█▊        | 209/1165 [00:02<00:12, 74.15it/s]
 19%|█▉        | 220/1165 [00:02<00:11, 80.19it/s]
 20%|██        | 234/1165 [00:02<00:09, 94.87it/s]
 21%|██        | 244/1165 [00:02<00:10, 88.09it/s]
 22%|██▏       | 254/1165 [00:02<00:12, 71.89it/s]
 23%|██▎       | 268/1165 [00:02<00:10, 87.04it/s]
 24%|██▍       | 281/1165 [00:03<00:09, 93.03it/s]
 25%|██▌       | 292/1165 [00:03<00:09, 96.61it/s]
 26%|██▌       | 303/1165 [00:03<00:09, 90.01it/s]
 28%|██▊       | 324/1165 [00:03<00:07, 118.13it/s]
 29%|██▉       | 340/1165 [00:03<00:06, 127.69it/s]
 30%|███       | 354/1165 [00:03<00:07, 104.92it/s]
 31%|███▏      | 366/1165 [00:04<00:10, 76.38it/s] 
 32%|███▏      | 377/1165 [00:04<00:09, 81.88it/s]
 33%|███▎      | 387/1165 [00:04<00:09, 79.50it/s]
 35%|███▍      | 404/1165 [00:04<00:07, 99.05it/s]
 36%|███▌      | 416/1165 [00:04<00:07, 100.76it/s]
 37%|███▋      | 428/1165 [00:04<00:08, 91.70it/s] 
 38%|███▊      | 438/1165 [00:04<00:08, 88.14it/s]
 39%|███▊      | 450/1165 [00:04<00:07, 94.92it/s]
 40%|███▉      | 461/1165 [00:05<00:08, 87.91it/s]
 41%|████      | 475/1165 [00:05<00:07, 97.71it/s]
 42%|████▏     | 486/1165 [00:05<00:07, 90.98it/s]
 43%|████▎     | 496/1165 [00:05<00:07, 92.22it/s]
 43%|████▎     | 506/1165 [00:05<00:07, 83.89it/s]
 45%|████▍     | 521/1165 [00:05<00:06, 99.24it/s]
 46%|████▌     | 532/1165 [00:05<00:08, 76.21it/s]
 46%|████▋     | 541/1165 [00:05<00:08, 71.42it/s]
 47%|████▋     | 549/1165 [00:06<00:08, 71.09it/s]
 48%|████▊     | 560/1165 [00:06<00:08, 75.28it/s]
 49%|████▉     | 568/1165 [00:06<00:09, 62.02it/s]
 50%|████▉     | 580/1165 [00:06<00:07, 74.30it/s]
 51%|█████     | 593/1165 [00:06<00:06, 86.86it/s]
 52%|█████▏    | 603/1165 [00:06<00:07, 74.31it/s]
 53%|█████▎    | 612/1165 [00:06<00:07, 70.86it/s]
 53%|█████▎    | 620/1165 [00:07<00:08, 66.32it/s]
 54%|█████▍    | 628/1165 [00:07<00:07, 67.22it/s]
 55%|█████▍    | 636/1165 [00:07<00:08, 65.49it/s]
 56%|█████▌    | 647/1165 [00:07<00:07, 70.01it/s]
 56%|█████▌    | 655/1165 [00:07<00:07, 67.49it/s]
 57%|█████▋    | 663/1165 [00:07<00:07, 69.12it/s]
 58%|█████▊    | 676/1165 [00:07<00:05, 82.86it/s]
 59%|█████▉    | 689/1165 [00:07<00:05, 94.49it/s]
 60%|██████    | 699/1165 [00:08<00:04, 93.52it/s]
 61%|██████    | 709/1165 [00:08<00:05, 85.12it/s]
 62%|██████▏   | 718/1165 [00:08<00:05, 79.42it/s]
 62%|██████▏   | 727/1165 [00:08<00:05, 77.17it/s]
 63%|██████▎   | 735/1165 [00:08<00:06, 70.79it/s]
 64%|██████▍   | 743/1165 [00:08<00:06, 66.25it/s]
 65%|██████▍   | 753/1165 [00:08<00:05, 74.15it/s]
 66%|██████▌   | 767/1165 [00:08<00:04, 90.67it/s]
 67%|██████▋   | 777/1165 [00:09<00:04, 86.06it/s]
 67%|██████▋   | 786/1165 [00:09<00:05, 69.10it/s]
 68%|██████▊   | 794/1165 [00:09<00:05, 63.83it/s]
 69%|██████▉   | 802/1165 [00:09<00:05, 67.20it/s]
 70%|██████▉   | 811/1165 [00:09<00:04, 71.99it/s]
 70%|███████   | 819/1165 [00:09<00:04, 73.82it/s]
 71%|███████   | 828/1165 [00:09<00:04, 76.05it/s]
 72%|███████▏  | 836/1165 [00:10<00:05, 64.66it/s]
 73%|███████▎  | 848/1165 [00:10<00:04, 77.73it/s]
 74%|███████▎  | 857/1165 [00:10<00:04, 68.31it/s]
 74%|███████▍  | 865/1165 [00:10<00:04, 63.66it/s]
 75%|███████▍  | 872/1165 [00:10<00:05, 55.24it/s]
 75%|███████▌  | 879/1165 [00:10<00:05, 56.62it/s]
 76%|███████▌  | 885/1165 [00:10<00:05, 49.65it/s]
 76%|███████▋  | 891/1165 [00:11<00:06, 42.48it/s]
 77%|███████▋  | 896/1165 [00:11<00:06, 43.98it/s]
 77%|███████▋  | 901/1165 [00:11<00:06, 39.32it/s]
 78%|███████▊  | 906/1165 [00:11<00:06, 40.87it/s]
 79%|███████▉  | 920/1165 [00:11<00:04, 58.51it/s]
 79%|███████▉  | 926/1165 [00:11<00:04, 52.78it/s]
 80%|████████  | 932/1165 [00:11<00:04, 52.13it/s]
 81%|████████  | 938/1165 [00:12<00:04, 45.78it/s]
 81%|████████  | 943/1165 [00:12<00:04, 44.48it/s]
 81%|████████▏ | 948/1165 [00:12<00:05, 43.01it/s]
 82%|████████▏ | 953/1165 [00:12<00:05, 38.45it/s]
 82%|████████▏ | 957/1165 [00:12<00:05, 37.33it/s]
 82%|████████▏ | 961/1165 [00:12<00:05, 37.27it/s]
 83%|████████▎ | 965/1165 [00:12<00:05, 37.26it/s]
 83%|████████▎ | 969/1165 [00:12<00:05, 33.19it/s]
 84%|████████▎ | 975/1165 [00:13<00:04, 39.74it/s]
 84%|████████▍ | 980/1165 [00:13<00:04, 40.05it/s]
 85%|████████▍ | 985/1165 [00:13<00:04, 37.65it/s]
 85%|████████▍ | 990/1165 [00:13<00:04, 38.64it/s]
 85%|████████▌ | 994/1165 [00:13<00:04, 37.96it/s]
 86%|████████▌ | 999/1165 [00:13<00:04, 38.75it/s]
 86%|████████▌ | 1003/1165 [00:13<00:04, 38.56it/s]
 87%|████████▋ | 1008/1165 [00:13<00:03, 40.32it/s]
 87%|████████▋ | 1013/1165 [00:14<00:03, 40.22it/s]
 87%|████████▋ | 1018/1165 [00:14<00:03, 42.16it/s]
 88%|████████▊ | 1023/1165 [00:14<00:03, 40.10it/s]
 88%|████████▊ | 1029/1165 [00:14<00:03, 45.13it/s]
 89%|████████▉ | 1036/1165 [00:14<00:02, 50.25it/s]
 89%|████████▉ | 1042/1165 [00:14<00:02, 45.34it/s]
 90%|████████▉ | 1047/1165 [00:14<00:02, 41.17it/s]
 90%|█████████ | 1052/1165 [00:14<00:02, 39.43it/s]
 91%|█████████ | 1057/1165 [00:15<00:02, 40.95it/s]
 91%|█████████▏| 1064/1165 [00:15<00:02, 45.63it/s]
 92%|█████████▏| 1071/1165 [00:15<00:01, 49.27it/s]
 92%|█████████▏| 1077/1165 [00:15<00:01, 48.06it/s]
 93%|█████████▎| 1082/1165 [00:15<00:01, 47.94it/s]
 93%|█████████▎| 1087/1165 [00:15<00:01, 47.09it/s]
 94%|█████████▍| 1095/1165 [00:15<00:01, 55.77it/s]
 95%|█████████▍| 1101/1165 [00:15<00:01, 51.01it/s]
 95%|█████████▌| 1110/1165 [00:16<00:00, 59.49it/s]
 96%|█████████▌| 1117/1165 [00:16<00:00, 54.37it/s]
 96%|█████████▋| 1123/1165 [00:16<00:00, 49.88it/s]
 97%|█████████▋| 1129/1165 [00:16<00:00, 47.47it/s]
 98%|█████████▊| 1137/1165 [00:16<00:00, 52.87it/s]
 98%|█████████▊| 1143/1165 [00:16<00:00, 45.61it/s]
 99%|█████████▊| 1149/1165 [00:16<00:00, 45.09it/s]
 99%|█████████▉| 1154/1165 [00:17<00:00, 43.19it/s]
 99%|█████████▉| 1159/1165 [00:17<00:00, 44.21it/s]
100%|█████████▉| 1164/1165 [00:17<00:00, 41.13it/s]
100%|██████████| 1165/1165 [00:17<00:00, 67.46it/s]


  0%|          | 0/156 [00:00<?, ?it/s]
 99%|█████████▊| 154/156 [00:00<00:00, 1533.81it/s]
100%|██████████| 156/156 [00:00<00:00, 1518.02it/s]

2026-09-29 17:45:24 - INFO - Removed 67 isolated nodes.

edges["tf"].value_counts().head(10)
tf BCL11A 30 MEF2C 23 KLF4 11 POU2AF1 11 IRF8 8 IRF4 4 YBX3 3 EBF1 3 TCF4 3 CREB5 2 Name: count, dtype: int64

The network consists of 115 links from 23 TFs to 84 genes, most of them from regulators of B cells, pDCs and monocytes such as BCL11A, MEF2C, POU2AF1, IRF8 and KLF4. link_tf_to_gene only keeps TFs whose own gene is part of the network, so regulators that are not among the highly variable genes, such as LEF1, TCF7, CEBPA and SPI1, are missing although their motifs are enriched in the previous chapter.

Comparison with prior knowledge

To see how much of the network is supported by known regulation, we compare its TF-gene links with CollecTRI Müller-Dott et al., 2023, the prior-knowledge network that we use in the RNA section.

collectri = ln.Artifact.get(key="conditions/gsea_pathway/collectri.parquet").load()
known = set(zip(collectri["source"], collectri["target"], strict=True))
edges["in_collectri"] = [
    edge in known for edge in zip(edges["tf"], edges["target"], strict=True)
]
edges[edges["in_collectri"]]
Loading...

Only 4 of the 115 links are in CollecTRI, although CollecTRI contains all 23 TFs. Prior knowledge is incomplete and not specific to cell types, so the other links are not necessarily wrong, but they remain hypotheses until they are tested, for example by perturbing the TF.

Quiz

Theory

Loading...

Networks from multiome data

Loading...

Contributors

We gratefully acknowledge the contributions of:

Authors

  • Ignacio Ibarra

Reviewers

  • Lukas Heumos

  • Anna Schaar

References
  1. Spitz, F., & Furlong, E. E. M. (2012). Transcription factors: from enhancer binding to developmental control. Nat. Rev. Genet., 13(9), 613–626.
  2. Bravo González-Blas, C., De Winter, S., Hulselmans, G., Hecker, N., Matetovici, I., Christiaens, V., Poovathingal, S., Wouters, J., Aibar, S., & Aerts, S. (2023). SCENIC+: single-cell multiomic inference of enhancers and gene regulatory networks. Nature Methods, 20(9), 1355–1367. 10.1038/s41592-023-01938-4
  3. Badia-i Mompel, P., Casals-Franch, R., Wessels, L., Müller-Dott, S., Trimbour, R., Yang, Y., Ramirez Flores, R. O., & Saez-Rodriguez, J. (2024). Comparison and evaluation of methods to infer gene regulatory networks from multimodal single-cell data. bioRxiv. 10.1101/2024.12.20.629764
  4. Weirauch, M. T., Yang, A., Albu, M., Cote, A. G., Montenegro-Montero, A., Drewe, P., Najafabadi, H. S., Lambert, S. A., Mann, I., Cook, K., Zheng, H., Goity, A., van Bakel, H., Lozano, J.-C., Galli, M., Lewsey, M. G., Huang, E., Mukherjee, T., Chen, X., … Hughes, T. R. (2014). Determination and Inference of Eukaryotic Transcription Factor Sequence Specificity. Cell, 158(6), 1431–1443. 10.1016/j.cell.2014.08.009
  5. Zhang, K., Zemke, N. R., Armand, E. J., & Ren, B. (2024). A fast, scalable and versatile tool for analysis of single-cell omics data. Nature Methods, 21(2), 217–227. 10.1038/s41592-023-02139-9
  6. Müller-Dott, S., Tsirvouli, E., Vazquez, M., Ramirez Flores, R. O., Badia-i Mompel, P., Fallegger, R., Türei, D., Lægreid, A., & Saez-Rodriguez, J. (2023). Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities. Nucleic Acids Research, 51(20), 10934–10949. 10.1093/nar/gkad841
Chromatin Accessibility
Peak calling and motif analysis
Spatial omics
Single-cell data resolved in space