Skip to content
PIASO

cytorete: RNA regulon inference

cytorete infers regulons — a transcription factor with its weighted target set — from RNA alone, on the PIASO stack: promoter windows from the PIASO-data TSS annotation, sequence from the 2bit genome, motifs scanned with piaso.pp.scan_motifs, and edges kept only when the TF’s and the target’s COSG cell-type specificity profiles agree (piaso.tl.cospecificity_trans).

Terminal window
pip install cytorete        # pulls piaso-tools>=1.2.2, cosg, cytome

1. Inputs — all fetchable

import piaso


jaspar = piaso.data.fetch_jaspar()                # JASPAR CORE vertebrates MEME
twobit = piaso.data.fetch_2bit("hg38")            # genome sequence
piaso.data.fetch_genome("hg38")                   # TSS BED among the references
ds = piaso.data.load_dataset("sea_ad_mtg_20k_cytome",
                                      return_type="cytome")
piaso.settings.set_figure_params(style="cell")   # one house style across every figure

Human cortex, 20,000 nuclei, Subclass annotations — streamed from the .cytome, never loaded whole.

2. One call

import cytorete


TFS = ["SOX9", "OLIG2", "SPI1", "IRF8", "MEF2C", "NEUROD2", "DLX2",
       "LHX6", "TCF7L2", "NFIB", "RORB", "CUX2", "FEZF2", "PAX6", "EGR1"]


cytorete.inferRegulon(ds, "hg38", "Subclass",
                      jaspar_path=jaspar, twobit_path=twobit,
                      tf_list=TFS)

What happened, from the run log:

[inferRegulon] 2553 genes with promoters (5396 intervals)
[cistrome]     20 PWMs (18 TFs) × 5396 promoter sequences  (method=motif_bg)
[cistrome]     M = 18 TFs × 2553 genes, 3419 edges (7.4% density)
[regulons]     15 global regulons (median 45 targets); per-cell-type regulons for 23 cell types
[regulonActivity]     scoring 15 regulons (layer=infog, pvalues=True)
[regulonSpecificity]  15 regulons × 24 cell types

Everything is written back onto the object: per-cell activity in X_regulon (+ X_regulon_pval), the regulon store under the regulon key. Omit tf_list to run every expressed TF with a motif.

3. Read the biology

cytorete.pl.regulonActivity(ds, groupby="Subclass")
Regulon activity per subclass

The controls fall where they should without any tuning: SPI1/IRF8 peak in microglia, OLIG2 in OPCs and oligodendrocytes, SOX9/PAX6 in astrocytes, LHX6 in the MGE-derived Pvalb/Sst interneurons, FEZF2 in L5/6 projection neurons, NEUROD2 broadly excitatory, NFIB in astrocytes and endothelium.

The clustered specificity view — run here with a larger TF list — sharpens the same picture into a near-diagonal: each lineage claims its regulons (OLIG1/2 → OPC, DLX2 → CGE interneurons, LHX6/NKX2-1 → MGE, CUX2 → L2/3 IT, BCL11B/FEZF2 → L5 ET and L5/6 NP, TBR1 → L6 IT Car3):

cytorete.pl.regulonSpecificityScatter  # ranked per cell type
cytorete.regulonSpecificity(ds, groupby="Subclass")
Regulon specificity clustermap

Also available: pl.regulonEmbedding (activity on the UMAP/spatial embedding) and pl.regulonNetwork (the TF→target graph).

3a. Which cells are which

Every regulon panel below is drawn on this UMAP, so it is worth seeing the cell types on it first:

piaso.pl.plotEmbedding(ds, color="Subclass", basis="X_umap")
SEA-AD subclasses on the UMAP

3b. Activity comes with a p-value

compute_pvalues=True is the default, so alongside X_regulon there is X_regulon_pval, from the same control-set null that PIASOscore uses:

import numpy as np
pvals = np.asarray(ds.embeddings["X_regulon_pval"])   # cells × regulons
names = md["names"]
j = names.index("CUX2")
(pvals[:, j] < 0.05).mean()          # fraction of cells with significant CUX2

Plot it beside the activity — the pair is more informative than either:

import numpy as np
P = np.asarray(ds.embeddings["X_regulon_pval"])
for tf in ["CUX2", "FEZF2", "LHX6", "IRF8"]:
    j = names.index(tf)
    ds.cells[f"act_{tf}"] = np.asarray(ds.embeddings["X_regulon"])[:, j]
    ds.cells[f"nlp_{tf}"] = -np.log10(np.clip(P[:, j], 1e-300, 1))
ds.flush()


piaso.pl.plotEmbedding(ds, color="act_CUX2", basis="X_umap", vmin_pct=5, vmax_pct=95)
piaso.pl.plotEmbedding(ds, color="nlp_CUX2", basis="X_umap", vmin_pct=5, vmax_pct=95)
Regulon activity and its significance

Read the two rows together. CUX2, FEZF2 and LHX6 reach −log10 p = 3.0 across their territories — and that is the ceiling, not a coincidence: with 1,000 control sets the smallest attainable p is 1/1001, so anything strongly significant piles up at exactly 3.0. The panels are saturated, which is the correct reading of “as significant as this test can report”.

IRF8 is the instructive one. Its activity is unmistakable in microglia, yet only 6% of cells clear p < 0.05 against 42–57% for the other three. Nothing is wrong: microglia are a small fraction of this dataset, the regulon is inactive nearly everywhere else, and a per-cell test asks whether this cell exceeds its own controls. A regulon can be highly specific and significant in few cells — specificity and per-cell significance are different questions, and the two rows let you see which one you are looking at.

Plot it like any other per-cell value:

for tf in ["CUX2", "FEZF2", "LHX6", "IRF8"]:
    j = names.index(tf)
    ds.cells[f"negp_{tf}"] = -np.log10(np.clip(pvals[:, j], 1e-300, 1))
ds.flush()


piaso.pl.plotEmbedding(ds, color="negp_CUX2", basis="X_umap", vmin=0, vmax=3)
Per-cell significance of four regulons

Significance is not a copy of activity. FEZF2 reaches p < 0.05 in 57% of cells and IRF8 in 6% — IRF8’s regulon is small and confined to microglia, so almost everywhere else it is indistinguishable from its control sets, which is the correct answer rather than a weak one. Reading the activity panel alone would suggest a graded signal across the section; the p-value says where that gradient is worth interpreting.

The p-value is permutation-based, so it floors at 1/(n_permutations+1) — -log10(p) saturates at 3.0 with the default 1,000 control sets, and values at the ceiling are “as significant as this test can report”, not “more significant than the others there.”

3c. The other three views

cytorete.pl has four plot functions; the heatmap above is one. The other three answer different questions on the same result.

Where is a regulon active? — activity on any embedding:

cytorete.pl.regulonEmbedding(ds, regulons=["CUX2", "FEZF2", "LHX6", "IRF8"],
                             basis="X_umap", key="X_regulon")
Regulon activity on the UMAP

Four regulons, four territories, none of them supplied to the method: CUX2 in upper-layer excitatory neurons, FEZF2 in deep layers, LHX6 in MGE-derived interneurons, IRF8 in microglia.

What is in a regulon? — the TF and its targets:

cytorete.pl.regulonNetwork(ds, tf="CUX2", max_targets=20)
CUX2 regulon network

Target names are drawn by default (label_targets=False for dense multi-TF layouts). CUX2’s here are synaptic and channel genes — GRIA3, KCNQ5, DLGAP1, CDH9 — which is the sanity check worth doing before believing a regulon.

Is the TF’s own expression driving it? — regulonSpecificityScatter plots regulon specificity against TF expression per cell type, which separates a TF that is specific because it is expressed there from one whose targets are coherent there.

What is not here yet

This release infers regulons from RNA. The multiome route — cistromes built from accessible peaks rather than promoter windows alone, using paired RNA + ATAC — is in development and will be released separately. Its entry points (inferGRN, inferTFActivity) exist in the package today and raise an ImportError naming what they need, so import cytorete behaves the same either way.

4. How the edges are decided

  1. Promoter cistrome — strand-aware windows (−1000/+500 around each TSS; alternative promoters kept separate) scanned against the motif DB with a per-motif background model, giving TF → gene motif support.
  2. Trans co-specificity — each motif-supported pair is scored by the cosine agreement of the TF’s and target’s COSG specificity profiles across your groupby cell types; positive-sign edges above the threshold survive, per cell type.
  3. Regulon assembly keeps TFs with at least min_targets surviving targets. NES pruning is available — cistrome_method="nes" switches the cistrome to an RcisTarget-style normalised enrichment score with nes_threshold — but it is not the default; motif_bg is.
  4. Activity is piaso.tl.score (PIASOscore), PIASO’s own gene-set scorer, not AUCell: for each cell it scores the regulon against size- and expression-matched control gene sets drawn from the same data, so a large regulon and a small one are comparable. Because the controls give a null, compute_pvalues=True (the default) also writes a per-cell p-value beside the score.

Works identically on an AnnData; on a cytome every stage streams. The motif engine itself — PWMs, backgrounds, p-value→score thresholds — is documented in the Motif analysis tutorial in this section.