Contribution maps: which bases drive accessibility¶

A contribution (attribution) map scores every base of a sequence by how much it moves a PEAgent readout. peagent_tool computes them by bottleneck in-silico mutagenesis: every position is mutated to all four nucleotides, the change in the shared 64-d embedding is measured once, and that single result is projected onto whichever cell group you ask for. Asking for one cell type and asking for all 23 cell types of a tissue therefore cost the same.

Two normalizations are available, and they answer different questions:

normalize	What it shows	Matches
"none"	Total drive on the target readout	PEAgent global browser track
"tissue"	Drive specific to this cell type relative to the other cell types of the same tissue	PEAgent cell-type browser tracks, and the TF-MoDISco input used in the paper

Everything below is real executed output.

In [1]:

from pathlib import Path
from urllib.request import urlretrieve

import numpy as np
import pandas as pd

import peagent_tool

print("peagent_tool", peagent_tool.__version__)

# The two example ACR windows used throughout this tutorial.
FASTA = Path("maize_example_acrs.fa")
if not FASTA.exists():
    urlretrieve(
        "https://raw.githubusercontent.com/YAOJ-bioin/peagent_tool/main/docs/assets/tutorials/maize_example_acrs.fa",
        FASTA,
    )
print(FASTA.read_text().splitlines()[0])

from pathlib import Path from urllib.request import urlretrieve import numpy as np import pandas as pd import peagent_tool print("peagent_tool", peagent_tool.__version__) # The two example ACR windows used throughout this tutorial. FASTA = Path("maize_example_acrs.fa") if not FASTA.exists(): urlretrieve( "https://raw.githubusercontent.com/YAOJ-bioin/peagent_tool/main/docs/assets/tutorials/maize_example_acrs.fa", FASTA, ) print(FASTA.read_text().splitlines()[0])

peagent_tool 0.1.0
>chr1_284657363_284657864 maize B73 NAM-5.0 | crownRoot_trichoblast cell-type-specific ACR | 1344 bp model window | window_start=284656941

1. Build the scorer¶

ContributionScorer loads the model once. Reuse it for every sequence.

In [2]:

from peagent_tool import ContributionScorer

scorer = ContributionScorer("maize")
print("tissues:", scorer.tissues)
print("cell types in crownRoot:", len(scorer.celltypes_in_tissue("crownRoot")))
scorer.celltypes_in_tissue("crownRoot")[:6]

from peagent_tool import ContributionScorer scorer = ContributionScorer("maize") print("tissues:", scorer.tissues) print("cell types in crownRoot:", len(scorer.celltypes_in_tissue("crownRoot"))) scorer.celltypes_in_tissue("crownRoot")[:6]

tissues: ['axillaryBud', 'crownRoot', 'ear', 'leaf', 'root', 'tassel']
cell types in crownRoot: 23

Out[2]:

['crownRoot_procambial_phloem_SE_precursors',
 'crownRoot_cortex',
 'crownRoot_vascular_sclerenchyma',
 'crownRoot_pith_parenchyma',
 'crownRoot_vascular_parenchyma',
 'crownRoot_trichoblast']

In [3]:

records = {}
for block in FASTA.read_text().split(">")[1:]:
    header, _, body = block.partition("\n")
    records[header.split()[0]] = "".join(body.split())

CTACR = "chr1_284657363_284657864"       # crownRoot_trichoblast cell-type-specific ACR
TARGET = "crownRoot_trichoblast"

# The ACR itself occupies the central 501 bp of the 1344 bp model window.
ACR_START, ACR_END = 423, 923
sequence = records[CTACR]
len(sequence)

records = {} for block in FASTA.read_text().split(">")[1:]: header, _, body = block.partition("\n") records[header.split()[0]] = "".join(body.split()) CTACR = "chr1_284657363_284657864" # crownRoot_trichoblast cell-type-specific ACR TARGET = "crownRoot_trichoblast" # The ACR itself occupies the central 501 bp of the 1344 bp model window. ACR_START, ACR_END = 423, 923 sequence = records[CTACR] len(sequence)

Out[3]:

1344

2. The global map (raw scores)¶

target="global" averages the readout over all 50,639 cells, which is the bulk-accessibility view. normalize="none" keeps the raw scores.

In [4]:

global_map = scorer.compute({CTACR: sequence}, target="global", normalize="none")
pd.Series(global_map.describe())

global_map = scorer.compute({CTACR: sequence}, target="global", normalize="none") pd.Series(global_map.describe())

Out[4]:

sequence_id                  chr1_284657363_284657864
species                                         maize
target                                         global
tissue                                           None
celltype                                         None
n_cells                                         50639
normalize                                        none
n_tissue_reference_groups                           0
normalized_length                                1344
normalization_action                        unchanged
dtype: object

In [5]:

from peagent_tool import plot_attribution

_ = plot_attribution(global_map, plot_start=ACR_START, plot_end=ACR_END,
                 title="Maize | global | raw contribution")

from peagent_tool import plot_attribution _ = plot_attribution(global_map, plot_start=ACR_START, plot_end=ACR_END, title="Maize | global | raw contribution")

3. One cell type, raw vs tissue-normalized¶

Now target a single tissue_celltype. The raw map is dominated by the same tissue-wide signal visible in the global track. Subtracting the mean over all 23 crownRoot cell types isolates what is specific to trichoblast — this is the version plotted as Attribution score in the PEAgent browser.

In [6]:

raw = scorer.compute({CTACR: sequence}, target=TARGET, normalize="none")
normalized = scorer.compute({CTACR: sequence}, target=TARGET, normalize="tissue")

print("groups averaged for the tissue reference:", len(normalized.tissue_reference))
print("first three:", normalized.tissue_reference[:3])
pd.DataFrame({
    "raw": pd.Series(raw.describe()),
    "tissue-normalized": pd.Series(normalized.describe()),
})

raw = scorer.compute({CTACR: sequence}, target=TARGET, normalize="none") normalized = scorer.compute({CTACR: sequence}, target=TARGET, normalize="tissue") print("groups averaged for the tissue reference:", len(normalized.tissue_reference)) print("first three:", normalized.tissue_reference[:3]) pd.DataFrame({ "raw": pd.Series(raw.describe()), "tissue-normalized": pd.Series(normalized.describe()), })

groups averaged for the tissue reference: 23
first three: ('crownRoot_procambial_phloem_SE_precursors', 'crownRoot_cortex', 'crownRoot_vascular_sclerenchyma')

Out[6]:

	raw	tissue-normalized
sequence_id	chr1_284657363_284657864	chr1_284657363_284657864
species	maize	maize
target	crownRoot_trichoblast	crownRoot_trichoblast
tissue	crownRoot	crownRoot
celltype	trichoblast	trichoblast
n_cells	370	370
normalize	none	tissue
n_tissue_reference_groups	0	23
normalized_length	1344	1344
normalization_action	unchanged	unchanged

In [7]:

from peagent_tool import plot_attribution_panels

_ = plot_attribution_panels(
    [raw, normalized],
    labels=["raw", "tissue-\nnormalized"],
    plot_start=ACR_START,
    plot_end=ACR_END,
    share_y=False,
    panel_height=1.5,
    title="Maize | %s" % TARGET,
)

from peagent_tool import plot_attribution_panels _ = plot_attribution_panels( [raw, normalized], labels=["raw", "tissue-\nnormalized"], plot_start=ACR_START, plot_end=ACR_END, share_y=False, panel_height=1.5, title="Maize | %s" % TARGET, )

The two versions agree where the tissue as a whole is driven, and diverge where trichoblast differs from its neighbours — here the normalized amplitude is about a third of the raw one:

In [8]:

window = slice(ACR_START - 1, ACR_END)
pd.Series({
    "pearson r (raw vs normalized)": round(float(np.corrcoef(
        raw.profile[window], normalized.profile[window])[0, 1]), 4),
    "max |raw|": round(float(np.abs(raw.profile[window]).max()), 4),
    "max |normalized|": round(float(np.abs(normalized.profile[window]).max()), 4),
})

window = slice(ACR_START - 1, ACR_END) pd.Series({ "pearson r (raw vs normalized)": round(float(np.corrcoef( raw.profile[window], normalized.profile[window])[0, 1]), 4), "max |raw|": round(float(np.abs(raw.profile[window]).max()), 4), "max |normalized|": round(float(np.abs(normalized.profile[window]).max()), 4), })

Out[8]:

pearson r (raw vs normalized)    0.5883
max |raw|                        0.3571
max |normalized|                 0.1174
dtype: float64

4. All cell types of a tissue in one pass¶

compute_many shares the single ISM pass across targets. Below, all 23 crownRoot cell types are scored and the strongest position of the normalized map is used to rank them — the cell type this ACR was called specific to should come out on top.

In [9]:

crownroot = scorer.celltypes_in_tissue("crownRoot")

all_maps = scorer.compute_many({CTACR: sequence}, targets=crownroot, normalize="tissue")
ranking = pd.DataFrame(
    [
        {
            "tissue_celltype": group,
            "n_cells": result.n_cells,
            "max_positive_contribution": float(result.profile[window].max()),
            "sum_positive_contribution": float(result.profile[window].clip(min=0).sum()),
        }
        for group, result in all_maps.items()
    ]
).sort_values("sum_positive_contribution", ascending=False).reset_index(drop=True)
ranking.head(8)

crownroot = scorer.celltypes_in_tissue("crownRoot") all_maps = scorer.compute_many({CTACR: sequence}, targets=crownroot, normalize="tissue") ranking = pd.DataFrame( [ { "tissue_celltype": group, "n_cells": result.n_cells, "max_positive_contribution": float(result.profile[window].max()), "sum_positive_contribution": float(result.profile[window].clip(min=0).sum()), } for group, result in all_maps.items() ] ).sort_values("sum_positive_contribution", ascending=False).reset_index(drop=True) ranking.head(8)

Out[9]:

	tissue_celltype	n_cells	max_positive_contribution	sum_positive_contribution
0	crownRoot_phloem_sieve_element_precursors	427	0.109990	9.347403
1	crownRoot_QC	590	0.110990	9.274614
2	crownRoot_columella/epidermal_initials	318	0.116833	8.919630
3	crownRoot_companion_cell	184	0.124925	8.607933
4	crownRoot_trichoblast	370	0.117407	8.427594
5	crownRoot_proximal_meristem	328	0.090682	7.402743
6	crownRoot_procambial_phloem_SE_precursors	471	0.087922	6.717281
7	crownRoot_cortex_endodermis_initials	407	0.078343	4.843545

Stacking the top cell types over the same window shows where they agree and where they differ. This is the same layout as the cell-type contribution panel in the PEAgent browser.

In [10]:

peak = int(np.argmax(normalized.profile[window])) + ACR_START
zoom_start, zoom_end = max(1, peak - 45), min(1344, peak + 45)
print("strongest cortex-specific position in the window:", peak)

top_groups = ranking["tissue_celltype"].head(6).tolist()
_ = plot_attribution_panels(
    [all_maps[group] for group in top_groups],
    labels=[group.replace("crownRoot_", "") for group in top_groups],
    plot_start=zoom_start,
    plot_end=zoom_end,
    panel_height=1.0,
    title="crownRoot cell types | tissue-normalized",
)

peak = int(np.argmax(normalized.profile[window])) + ACR_START zoom_start, zoom_end = max(1, peak - 45), min(1344, peak + 45) print("strongest cortex-specific position in the window:", peak) top_groups = ranking["tissue_celltype"].head(6).tolist() _ = plot_attribution_panels( [all_maps[group] for group in top_groups], labels=[group.replace("crownRoot_", "") for group in top_groups], plot_start=zoom_start, plot_end=zoom_end, panel_height=1.0, title="crownRoot cell types | tissue-normalized", )

strongest cortex-specific position in the window: 684

5. Saving publication files¶

Passing out_prefix writes a vector PDF and a 600 dpi PNG with Illustrator-compatible fonts.

In [11]:

paths = plot_attribution(
    normalized,
    out_prefix="crownRoot_trichoblast_contribution",
    plot_start=zoom_start,
    plot_end=zoom_end,
)
{key: str(value) for key, value in paths.items()}

paths = plot_attribution( normalized, out_prefix="crownRoot_trichoblast_contribution", plot_start=zoom_start, plot_end=zoom_end, ) {key: str(value) for key, value in paths.items()}

Out[11]:

{'pdf': 'crownRoot_trichoblast_contribution.pdf',
 'png': 'crownRoot_trichoblast_contribution.png',
 'font': 'Arial Unicode MS'}

6. Tidy output for downstream analysis¶

to_dataframe gives a per-position table carrying both score versions, and position_offset converts model-window coordinates into genome coordinates.

In [12]:

WINDOW_START = 284656941   # from the FASTA header
table = normalized.to_dataframe(position_offset=WINDOW_START)
table.iloc[peak - 4 : peak + 4]

WINDOW_START = 284656941 # from the FASTA header table = normalized.to_dataframe(position_offset=WINDOW_START) table.iloc[peak - 4 : peak + 4]

Out[12]:

	position	base	contribution	contribution_raw
680	284657622	T	0.061253	0.123613
681	284657623	A	0.027607	0.061933
682	284657624	T	0.079295	-0.012182
683	284657625	C	0.117407	-0.061837
684	284657626	G	-0.005414	0.029296
685	284657627	C	0.046760	0.080757
686	284657628	C	-0.039626	0.048050
687	284657629	C	0.006860	0.059300

The full A/C/G/T matrices are available too: attribution is the reference-projected map drawn above, and hypothetical keeps the score of every possible base at every position, which is what motif-discovery tools such as TF-MoDISco consume.

In [13]:

print("attribution ", normalized.attribution.shape)
print("hypothetical", normalized.hypothetical.shape)
normalized.hypothetical_dataframe().iloc[peak - 2 : peak + 2]

print("attribution ", normalized.attribution.shape) print("hypothetical", normalized.hypothetical.shape) normalized.hypothetical_dataframe().iloc[peak - 2 : peak + 2]

attribution  (1344, 4)
hypothetical (1344, 4)

Out[13]:

	A	C	G	T
682	-0.012088	0.007600	-0.074807	0.079295
683	-0.082004	0.117407	-0.122879	0.087475
684	-0.054967	0.022858	-0.005414	0.037523
685	-0.051820	0.046760	-0.039711	0.044771

7. Synthetic sequences¶

Contribution maps do not need a genome. Implant a motif into a random background and read off what the model does with it — the basis for high-throughput screening of designed regulatory sequences.

In [14]:

rng = np.random.default_rng(2026)
background = "".join(rng.choice(list("ACGT"), size=1344, p=[0.29, 0.21, 0.21, 0.29]))
motif = "AAAGCAAAGC"
designed = background[:667] + motif + background[677:]

designed_map = scorer.compute({"designed": designed}, target=TARGET, normalize="none")
background_map = scorer.compute({"background": background}, target=TARGET, normalize="none")
pd.Series({
    "max |contribution| background": round(float(np.abs(background_map.profile[642:702]).max()), 4),
    "max |contribution| designed": round(float(np.abs(designed_map.profile[642:702]).max()), 4),
})

rng = np.random.default_rng(2026) background = "".join(rng.choice(list("ACGT"), size=1344, p=[0.29, 0.21, 0.21, 0.29])) motif = "AAAGCAAAGC" designed = background[:667] + motif + background[677:] designed_map = scorer.compute({"designed": designed}, target=TARGET, normalize="none") background_map = scorer.compute({"background": background}, target=TARGET, normalize="none") pd.Series({ "max |contribution| background": round(float(np.abs(background_map.profile[642:702]).max()), 4), "max |contribution| designed": round(float(np.abs(designed_map.profile[642:702]).max()), 4), })

Out[14]:

max |contribution| background    0.1675
max |contribution| designed      0.2301
dtype: float64

In [15]:

_ = plot_attribution_panels(
    [background_map, designed_map],
    labels=["background", "+ motif"],
    plot_start=642,
    plot_end=702,
    panel_height=1.3,
    title="Implanted %s | %s" % (motif, TARGET),
)

_ = plot_attribution_panels( [background_map, designed_map], labels=["background", "+ motif"], plot_start=642, plot_end=702, panel_height=1.3, title="Implanted %s | %s" % (motif, TARGET), )

Consistency with the PEAgent browser¶

The tissue-normalized scores computed here reproduce the published PEAgent attribution tracks; see Consistency check for the numbers.

Next¶

• Variant effects — score REF/ALT alleles of SNVs.
• In-silico synergy — test whether two motifs act more than additively.