Motif analysis: sequence, PWMs and what a hit is worth

Scanning is the easy half. piaso.data fetches the genome and the motif databases, piaso.pp.scan_motifs finds the hits, and the Rust backend does it 91× faster than numpy. The half that decides whether the answer means anything is what you compare against — and this page contains a result that looks completely convincing and is entirely an artefact.

Everything here is in pip install piaso-tools. Reading sequence from a .2bit needs one optional package:

# pip install py2bit
import numpy as np
import pandas as pd
import piaso
from scipy.stats import fisher_exact
from statsmodels.stats.multitest import multipletests

piaso.settings.set_figure_params(style="cell")

1. The genome

piaso.data keeps a small cache of genome files and tells you what it has:

piaso.data.list_available_genomes()

['hg38', 'mm10']

files = piaso.data.resolve_genome_files("hg38")
sorted(files)

['chrom_sizes', 'ctcf', 'gene_boundary', 'gtf', 'promoter', 'tss_bed']

tss_bed is one row per transcript — 146,857 of them — which is what promoter windows are built from.

tss = pd.read_csv(files["tss_bed"], sep="\t", header=None,
                  names=["chrom", "start", "end", "gene", "score",
                         "strand", "biotype", "gene_id", "tx_id"])
tss.shape

(146857, 9)

2. Promoter windows

TSS − 1000 to TSS + 500, on the gene’s strand. A minus-strand gene’s upstream is to the right in genome coordinates; getting this backwards silently scans the wrong 1.5 kb and produces a clean-looking wrong answer.

UP, DOWN = 1000, 500

def promoters(genes):
    """One window per gene, from its first listed TSS."""
    sub = tss[tss["gene"].isin(genes)].drop_duplicates("gene")
    rows = []
    for r in sub.itertuples():
        lo, hi = ((r.start - UP, r.start + DOWN) if r.strand == "+"
                  else (r.start - DOWN, r.start + UP))
        rows.append((r.gene, r.chrom, max(0, lo), hi, r.strand))
    return pd.DataFrame(rows, columns=["gene", "chrom", "start", "end", "strand"])

3. Gene sets: human brain cell classes

Marker sets from an atlas beat a hand-picked list — they are bigger, and they were not chosen by someone who already had a hypothesis about which motifs would win. SEA-AD middle temporal gyrus gives 24 subclasses, 50 markers each:

markers_df, marker_db = piaso.tl.getMarkers(study="SEAAD2024_MTG_Subclass",
                                            as_dict=True)
len(marker_db), sorted(marker_db)[:6]

(24, ['Astrocyte', 'Chandelier', 'Endothelial', 'L2/3 IT', 'L4 IT', 'L5 ET'])

Group them into classes, so each class has enough promoters to test:

CLASSES = {
    "Glia": ["Astrocyte", "Oligodendrocyte", "OPC", "Microglia-PVM"],
    "Excitatory": ["L2/3 IT", "L4 IT", "L5 IT", "L6 IT", "L5 ET", "L6 CT",
                   "L5/6 NP", "L6b", "L6 IT Car3"],
    "Inhibitory": ["Pvalb", "Sst", "Vip", "Lamp5", "Sncg", "Chandelier",
                   "Pax6", "Lamp5 Lhx6"],
    "Vascular": ["Endothelial", "VLMC"],
}
sets = {name: sorted({g for t in types for g in marker_db[t]})
        for name, types in CLASSES.items()}
{k: len(v) for k, v in sets.items()}

{'Glia': 200, 'Excitatory': 432, 'Inhibitory': 349, 'Vascular': 100}

4. Sequence

twobit = piaso.data.fetch_2bit("hg38")     # ~800 MB, cached after the first call

def seqs_for(df):
    return piaso.data.extract_sequences(
        twobit, [(r.chrom, r.start, r.end, r.strand) for r in df.itertuples()])

fg_seq = {k: seqs_for(promoters(v)) for k, v in sets.items()}
{k: len(v) for k, v in fg_seq.items()}

{'Glia': 200, 'Excitatory': 432, 'Inhibitory': 348, 'Vascular': 99}

extract_sequences takes the strand in the interval and reverse-complements for you, so every sequence reads 5′→3′ along the gene.

5. Motifs

pwms = piaso.data.load_jaspar_meme(piaso.data.fetch_jaspar())
p = pwms[0]
p.motif_id, p.tf_name, p.width, p.probs.shape

('MA0004.1', 'Arnt', 6, (4, 6))

Note the shape: probs is (4, width) — bases down the rows, positions across. Broadcasting a length-4 background against it needs bg[:, None], not bg[None, :].

build_tf_motif_map groups motifs by the TF that binds them and, given a gene universe, drops TFs absent from your data:

tf2pwm = piaso.data.build_tf_motif_map(pwms, tf_list=None,
                                       gene_universe=list(tss["gene"].unique()))
motifs = list({p.motif_id: p for ps in tf2pwm.values() for p in ps}.values())
len(tf2pwm), len(motifs)

(740, 795)

tf_list=None uses the TF names the motif database itself carries. Pass a curated list where you have one — piaso.data.load_tf_list reads an AnimalTFDB table, which piaso.data.fetch_animaltfdb_tf_list("human") downloads once.

6. Why the threshold is per motif

A score threshold is meaningless across motifs of different lengths, because a longer motif accumulates more score. pvalue_to_threshold converts a p-value into the score cutoff for that PWM against that background:

bg_freq = piaso.pp.estimate_background(fg_seq["Glia"])

for tf in ("NEUROD2", "SOX2", "CTCF"):
    p = tf2pwm[tf][0]
    pssm = np.log2((p.probs + 0.01) / bg_freq[:, None])
    thr = piaso.pp.pvalue_to_threshold(pssm, bg_freq, pvalue=1e-4)
    print(f"{tf:8s} {p.motif_id:9s} width {p.width:2d}  threshold {thr:.2f}")

CTCF is the longest and most informative motif of the three, and gets the lowest threshold. A fixed score cutoff would have called CTCF sites almost everywhere while missing NEUROD2 entirely. This is why scan_motifs takes pvalue=, not a score.

7. Scanning

hits = piaso.pp.scan_motifs(motifs, seqs, background=bg_freq, pvalue=1e-4)
{k: getattr(v, "shape", len(v)) for k, v in sorted(hits.items())}

{'best_score': (795, 1200), 'hit_count': (795, 1200),
 'motif_ids': 795, 'tf_names': 795}

Two motif × sequence matrices: the best score anywhere in each sequence, and how many positions passed. hit_count > 0 is “this promoter has a site”.

backend="auto" uses Rust when the extension is present. On 706 motifs × 1,212 promoters × 1.5 kb:

backend	time
numpy	113.6 s
rust	1.2 s

91×. piaso.pp.rust_ext_available() says whether you have it; the numpy path is the same contract, so nothing changes but the wait.

8. A convincing result that is not real

Test each motif for enrichment in the 200 glial marker promoters against 1,000 protein-coding promoters sampled at random.

def enrich(fg, bg):
    seqs = fg + bg
    hits = piaso.pp.scan_motifs(motifs, seqs,
                                background=piaso.pp.estimate_background(seqs),
                                pvalue=1e-4)
    present = np.asarray(hits["hit_count"]) > 0
    n, rows = len(fg), []
    for i, motif in enumerate(hits["motif_ids"]):
        a, b = int(present[i, :n].sum()), int(present[i, n:].sum())
        if a + b == 0:
            continue
        odds, pv = fisher_exact([[a, n - a], [b, len(bg) - b]],
                                alternative="greater")
        rows.append(dict(motif=motif, tf=hits["tf_names"][i],
                         fg_pct=100 * a / n, bg_pct=100 * b / len(bg),
                         odds=odds, pval=pv))
    out = pd.DataFrame(rows).sort_values("pval")
    out["qval"] = multipletests(out["pval"], method="fdr_bh")[1]
    return out

Glia vs random promoters: 4 of 789 at q<0.05

   motif      tf  fg_pct  bg_pct  odds     pval   qval
MA0087.3    SOX5    37.5    23.7  1.93 5.59e-05 0.0388
MA0041.3   FOXD3    75.5    61.9  1.90 0.000126 0.0388
MA0077.2    SOX9    31.0    19.1  1.90 0.000193 0.0388
MA1108.3    MXI1    60.0    46.0  1.76 0.000197 0.0388

Stop and read that. SOX9 is an astrocyte transcription factor. SOX5 is expressed across the glial lineage. FOXD3 is a neural-crest/glial-lineage factor. Three of the four top hits are the genes a reviewer would nod at. Nothing about this table looks wrong.

It is composition. Check it:

def gc(seqs):
    return np.array([(s.count("G") + s.count("C")) / len(s) for s in seqs])

gc(fg_seq["Glia"]).mean(), gc(pool_seq).mean()

(0.499, 0.528)

The glial promoters are 2.9 GC points poorer than the genomic pool. SOX motifs are AT-rich. An AT-rich motif finds more sites in AT-poorer sequence for no biological reason whatsoever — and it happens to name the right genes, because glial marker genes and SOX motifs are AT-rich for related but non-causal reasons.

Sample the background to match the foreground’s GC deciles instead:

edges = np.quantile(gc(fg), np.linspace(0, 1, 11))
edges[0], edges[-1] = 0.0, 1.0
fg_bin = np.digitize(gc(fg), edges[1:-1])
pool_bin = np.digitize(gc(pool_seq), edges[1:-1])

rng = np.random.default_rng(0)
picked = []
for b in range(11):
    want, have = int(5 * (fg_bin == b).sum()), np.flatnonzero(pool_bin == b)
    if want and have.size:
        picked.append(rng.choice(have, size=min(want, have.size), replace=False))
matched = np.concatenate(picked)

Run all four classes that way:

class	promoters	GC (fg / bg)	significant at q<0.05
Glia	200	0.499 / 0.499	0 of 790
Excitatory	432	0.466 / 0.467	0 of 785
Inhibitory	348	0.505 / 0.505	0 of 786
Vascular	99	0.525 / 0.526	0 of 787

Nothing survives, in any class. SOX9 and SOX5 fall out entirely. The four “discoveries” in the previous table were GC content wearing the names of plausible transcription factors.

The near-misses stay biologically sensible, which is worth noticing but not worth reporting as a finding: NEUROD1 and NEUROD2 are the 7th and 8th ranked motifs for excitatory neurons (q = 0.89), and LHX6 — the MGE interneuron factor — is 7th for inhibitory (q = 0.52). The signal is in the right direction and far too weak to call.

9. The contrast that does work

Testing marker promoters against any promoter asks “is being a marker gene associated with this motif”, which mixes the class with the fact of being a marker at all. The sharper question compares classes to each other: glial markers against neuronal markers, both sides marker promoters of comparable specificity, GC-matched.

neuron = sorted({g for t in CLASSES["Excitatory"] + CLASSES["Inhibitory"]
                 for g in marker_db[t]})

Glia: 200 promoters, GC 0.499
Neuron: 773 promoters, GC 0.484

Neuron > Glia: 773 vs 200 promoters, 2 of 795 at q<0.05

   motif    tf  fg_pct  bg_pct  odds     pval   qval
MA0613.1 FOXG1    87.2    74.5  2.33 1.86e-05  0.012
MA0842.3   NRL    92.4      82  2.66 3.01e-05  0.012
MA1106.2 HIF1A    48.1    34.5  1.76 0.000342 0.0756
MA0108.3   TBP    88.7      79  2.10 0.000392 0.0756
...
MA0885.3  DLX2    83.4      73.5 1.82  0.00123 0.0886

FOXG1, q = 0.012. FOXG1 is the forebrain neuronal transcription factor — it specifies telencephalic identity, and this is human middle temporal gyrus. DLX2, the GABAergic factor, sits just below the line at q = 0.089. Both are mechanistically right, and both survived the control that killed SOX9.

The reverse contrast, glia over neurons, gives nothing at q<0.05; MITF and NR2F1 lead at q = 1.

What to take from this

• Scanning is solved. 795 motifs across 1,200 promoters in about a second. It is not where the difficulty is.
• A motif hit is not a binding site. These promoters were scanned whether or not they are open in any cell. Sequence says where a TF could bind.
• Composition-matched backgrounds are not optional, and the unmatched result will not look wrong. It named SOX9 for glia. Match on GC and it disappears.
• The comparison matters more than the test. Four classes against a genomic background: zero. Two classes against each other: FOXG1, correctly.
• Even so, one motif out of 795 is a thin harvest from ~1,000 promoters. Promoter sequence alone is weak evidence for cell-type-specific regulation, because promoters are largely shared between cell types.

Which is the argument for using accessibility. Restricting the scan to peaks open in a given cell type, and requiring the TF’s own expression to track its targets’, is what turns a motif hit into a candidate regulatory edge — that is what cytorete is for, and it builds on exactly the functions above.

Regulatory regions beyond promoters

piaso.data also carries the SCREEN cCRE registry, so the same scan can run on candidate enhancers instead:

ccres = piaso.data.load_screen_ccres(piaso.data.fetch_screen("hg38"),
                                     classes=("PLS", "pELS", "dELS"))
sum(len(v["starts"]) for v in ccres.values())

row = promoters(["SOX9"]).iloc[0]
near = piaso.data.ccres_near_tss(ccres, row.chrom, int(row.start + UP), 100_000)
len(near)

Feed their sequences to extract_sequences and scan_motifs exactly as above — and match the background on GC, which matters more for enhancers than for promoters, not less.