In-silico synergy: do two motifs cooperate?¶

Two transcription-factor motifs may simply add up, or the pair may do more together than apart. peagent_tool tests this by in-silico marginalization: motif A alone, motif B alone and both together are implanted into a panel of background sequences over a grid of spacings and relative orientations, and the joint effect is compared with the additive expectation.

The readout is the pre-sigmoid pseudobulk logit of one cell group, f_g(x) = z(x) · w_g + b_g, and

synergy = Δjoint − (ΔA + ΔB)

is tested per arrangement with a one-sided paired Wilcoxon signed-rank test over backgrounds, then Benjamini-Hochberg corrected across all arrangements. This is the method of 05_synergy in the PEAgent analysis pipeline.

In [1]:

import numpy as np
import pandas as pd

import peagent_tool
from peagent_tool.synergy import MotifSynergyTester, random_backgrounds

print("peagent_tool", peagent_tool.__version__)

import numpy as np import pandas as pd import peagent_tool from peagent_tool.synergy import MotifSynergyTester, random_backgrounds print("peagent_tool", peagent_tool.__version__)

peagent_tool 0.1.0

1. Backgrounds and the tester¶

Backgrounds carry the sequence context against which the implanted motifs are measured. random_backgrounds draws GC-matched random sequence; shuffled_backgrounds dinucleotide-shuffles a real sequence when you want to preserve local composition. More backgrounds means a better-powered paired test.

In [2]:

backgrounds = random_backgrounds(n=48, gc_content=0.44, seed=2026)
tester = MotifSynergyTester("maize")

TARGET = "crownRoot_cortex"
weights, bias = tester.target_readout(TARGET)
baseline = tester.group_logit(backgrounds, weights, bias)
pd.Series({
    "backgrounds": len(backgrounds),
    "bottleneck size": len(weights),
    "baseline logit mean": round(float(baseline.mean()), 4),
    "baseline logit sd": round(float(baseline.std()), 4),
})

backgrounds = random_backgrounds(n=48, gc_content=0.44, seed=2026) tester = MotifSynergyTester("maize") TARGET = "crownRoot_cortex" weights, bias = tester.target_readout(TARGET) baseline = tester.group_logit(backgrounds, weights, bias) pd.Series({ "backgrounds": len(backgrounds), "bottleneck size": len(weights), "baseline logit mean": round(float(baseline.mean()), 4), "baseline logit sd": round(float(baseline.std()), 4), })

Out[2]:

backgrounds            48.0000
bottleneck size        64.0000
baseline logit mean    -4.1404
baseline logit sd       0.6503
dtype: float64

2. One pair, one arrangement grid¶

test_pair returns one row per (orientation, spacing) arrangement. Spacings that would make the two motifs overlap are skipped. The motifs below are a Dof-family core and a bZIP/G-box-like element, two families in the PEAgent contribution lexicon.

In [3]:

MOTIF_A = "AAAGCAAAG"    # Dof-family core
MOTIF_B = "TGACGTCA"     # bZIP / ACGT-core element

arrangements = tester.test_pair(
    MOTIF_A,
    MOTIF_B,
    target=TARGET,
    backgrounds=backgrounds,
    spacings=range(0, 41),
    name_a="Dof",
    name_b="bZIP",
)
print(arrangements.shape)
arrangements.sort_values("mean_synergy", ascending=False).head(8)[
    ["motif_a", "motif_b", "orientation", "spacing", "mean_delta_a", "mean_delta_b",
     "mean_delta_joint", "mean_synergy", "p"]
]

MOTIF_A = "AAAGCAAAG" # Dof-family core MOTIF_B = "TGACGTCA" # bZIP / ACGT-core element arrangements = tester.test_pair( MOTIF_A, MOTIF_B, target=TARGET, backgrounds=backgrounds, spacings=range(0, 41), name_a="Dof", name_b="bZIP", ) print(arrangements.shape) arrangements.sort_values("mean_synergy", ascending=False).head(8)[ ["motif_a", "motif_b", "orientation", "spacing", "mean_delta_a", "mean_delta_b", "mean_delta_joint", "mean_synergy", "p"] ]

(128, 11)

Out[3]:

	motif_a	motif_b	orientation	spacing	mean_delta_a	mean_delta_b	mean_delta_joint	mean_synergy	p
41	Dof	bZIP	+-	18	-0.01445	0.162772	0.351727	0.203405	0.001700
9	Dof	bZIP	++	18	-0.01445	0.162772	0.351727	0.203405	0.001700
8	Dof	bZIP	++	17	-0.01445	0.162772	0.324166	0.175844	0.005446
40	Dof	bZIP	+-	17	-0.01445	0.162772	0.324166	0.175844	0.005446
10	Dof	bZIP	++	19	-0.01445	0.162772	0.323224	0.174902	0.003156
42	Dof	bZIP	+-	19	-0.01445	0.162772	0.323224	0.174902	0.003156
43	Dof	bZIP	+-	20	-0.01445	0.162772	0.313056	0.164734	0.010461
11	Dof	bZIP	++	20	-0.01445	0.162772	0.313056	0.164734	0.010461

Each motif on its own has an effect; the question is whether the pair beats the sum. Positive mean_synergy with a small p-value means it does.

!!! note "Palindromic motifs" TGACGTCA is its own reverse complement, so the ++ and +- rows above are identical by construction, as are -+ and --. That is expected, not a bug — reverse-complementing a palindrome changes nothing. It also means a palindromic motif effectively halves the orientation grid.

In [4]:

pd.Series({
    "mean ΔA (Dof alone)": round(float(arrangements["mean_delta_a"].mean()), 4),
    "mean ΔB (bZIP alone)": round(float(arrangements["mean_delta_b"].mean()), 4),
    "best Δjoint": round(float(arrangements["mean_delta_joint"].max()), 4),
    "best synergy": round(float(arrangements["mean_synergy"].max()), 4),
    "arrangements tested": len(arrangements),
})

pd.Series({ "mean ΔA (Dof alone)": round(float(arrangements["mean_delta_a"].mean()), 4), "mean ΔB (bZIP alone)": round(float(arrangements["mean_delta_b"].mean()), 4), "best Δjoint": round(float(arrangements["mean_delta_joint"].max()), 4), "best synergy": round(float(arrangements["mean_synergy"].max()), 4), "arrangements tested": len(arrangements), })

Out[4]:

mean ΔA (Dof alone)      -0.0053
mean ΔB (bZIP alone)      0.1628
best Δjoint               0.3517
best synergy              0.2034
arrangements tested     128.0000
dtype: float64

3. Is the effect spacing-dependent?¶

Hard syntax — cooperation only at particular spacings — shows up as structure in the spacing profile. A flat profile means the two motifs merely co-occur.

In [5]:

import matplotlib.pyplot as plt
from matplotlib import font_manager

candidates = ["Arial", "Arial Unicode MS", "Liberation Sans", "Nimbus Sans", "DejaVu Sans"]
available = {font.name for font in font_manager.fontManager.ttflist}
plt.rcParams.update({
    "pdf.fonttype": 42, "ps.fonttype": 42, "font.family": "sans-serif",
    "font.sans-serif": [next((f for f in candidates if f in available), "DejaVu Sans")],
    "axes.labelsize": 14, "xtick.labelsize": 12, "ytick.labelsize": 12, "legend.fontsize": 11,
})

fig, ax = plt.subplots(figsize=(8, 4))
colors = {"++": "#D7191C", "+-": "#FDAE61", "-+": "#2C7BB6", "--": "#1A9850"}
for orientation, group in arrangements.groupby("orientation"):
    group = group.sort_values("spacing")
    ax.plot(group["spacing"], group["mean_synergy"], marker="o", markersize=3.5,
            color=colors.get(orientation, "black"), label=orientation)
ax.axhline(0.0, color="black", linewidth=1.0, linestyle="--")
ax.set_xlabel("Spacing between motif centres (bp)")
ax.set_ylabel("Mean synergy (logit)")
ax.grid(False)
for side in ("top", "right", "bottom", "left"):
    ax.spines[side].set_color("black")
    ax.spines[side].set_linewidth(1.2)
ax.legend(title="orientation", frameon=False, ncol=4)
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt from matplotlib import font_manager candidates = ["Arial", "Arial Unicode MS", "Liberation Sans", "Nimbus Sans", "DejaVu Sans"] available = {font.name for font in font_manager.fontManager.ttflist} plt.rcParams.update({ "pdf.fonttype": 42, "ps.fonttype": 42, "font.family": "sans-serif", "font.sans-serif": [next((f for f in candidates if f in available), "DejaVu Sans")], "axes.labelsize": 14, "xtick.labelsize": 12, "ytick.labelsize": 12, "legend.fontsize": 11, }) fig, ax = plt.subplots(figsize=(8, 4)) colors = {"++": "#D7191C", "+-": "#FDAE61", "-+": "#2C7BB6", "--": "#1A9850"} for orientation, group in arrangements.groupby("orientation"): group = group.sort_values("spacing") ax.plot(group["spacing"], group["mean_synergy"], marker="o", markersize=3.5, color=colors.get(orientation, "black"), label=orientation) ax.axhline(0.0, color="black", linewidth=1.0, linestyle="--") ax.set_xlabel("Spacing between motif centres (bp)") ax.set_ylabel("Mean synergy (logit)") ax.grid(False) for side in ("top", "right", "bottom", "left"): ax.spines[side].set_color("black") ax.spines[side].set_linewidth(1.2) ax.legend(title="orientation", frameon=False, ncol=4) fig.tight_layout() plt.show()

4. Several pairs with one correction¶

test_pairs runs a list of pairs and applies a single Benjamini-Hochberg correction across every arrangement of every pair, which is the correct unit of multiple testing when the spacing grid is searched. Homotypic pairs (motif_a == motif_b) automatically use the three distinguishable orientations.

In [6]:

pairs = [
    {"motif_a": "AAAGCAAAG", "motif_b": "TGACGTCA", "name_a": "Dof", "name_b": "bZIP"},
    {"motif_a": "AAAGCAAAG", "motif_b": "AAAGCAAAG", "name_a": "Dof", "name_b": "Dof"},
    {"motif_a": "TGACGTCA", "motif_b": "CCAAATTTGG", "name_a": "bZIP", "name_b": "CArG-like"},
]
result = tester.test_pairs(
    pairs,
    target=TARGET,
    backgrounds=backgrounds,
    spacings=range(0, 41, 2),
    alpha=0.01,
)
print("arrangements:", len(result.arrangements), " significant at padj<0.01:", len(result.significant(0.01)))
result.best[["motif_a", "motif_b", "orientation", "spacing", "mean_delta_joint", "mean_synergy", "p", "padj"]]

pairs = [ {"motif_a": "AAAGCAAAG", "motif_b": "TGACGTCA", "name_a": "Dof", "name_b": "bZIP"}, {"motif_a": "AAAGCAAAG", "motif_b": "AAAGCAAAG", "name_a": "Dof", "name_b": "Dof"}, {"motif_a": "TGACGTCA", "motif_b": "CCAAATTTGG", "name_a": "bZIP", "name_b": "CArG-like"}, ] result = tester.test_pairs( pairs, target=TARGET, backgrounds=backgrounds, spacings=range(0, 41, 2), alpha=0.01, ) print("arrangements:", len(result.arrangements), " significant at padj<0.01:", len(result.significant(0.01))) result.best[["motif_a", "motif_b", "orientation", "spacing", "mean_delta_joint", "mean_synergy", "p", "padj"]]

arrangements: 176  significant at padj<0.01: 0

Out[6]:

	motif_a	motif_b	orientation	spacing	mean_delta_joint	mean_synergy	p	padj
0	Dof	bZIP	+-	18	0.351727	0.203405	0.001700	0.149556
1	bZIP	CArG-like	+-	16	0.282628	0.138828	0.023845	0.225972
2	Dof	Dof	+-	18	0.117762	0.128457	0.022121	0.225972

!!! note "Degenerate motif strings" Motifs are implanted verbatim. Use a concrete consensus (TGACGTCA), not an IUPAC string with W/N, unless you intend those literal characters to be written into the background — N is encoded as a uniform 0.25 vector by the model and will damp the effect.

5. What to report¶

• mean_delta_joint — the size of the combined effect on the group logit.
• mean_synergy — how much of it exceeds the additive expectation.
• padj — significance after correction across all arrangements.
• result.best — the strongest arrangement per pair, one row per pair.

Next¶

• Contribution maps — find candidate motifs to test.
• Variant effects — score natural alleles.