Predicting cell-type chromatin accessibility from a sequence¶

This tutorial takes a DNA sequence and returns, for every cell type in a PEAgent atlas, the model's predicted probability that the sequence is accessible in that cell type.

Everything below is executed against the maize atlas model (50,639 nuclei, 94 tissue_celltype groups, B73 NAM-5.0); the same three lines work for soybean and rice. All output on this page is the real output of the code above it.

What you need:

• peagent_tool installed with the prediction extra (see Installation)
• the maize model weights and cell-group metadata in place (see Models and data)

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. Check that the model assets are in place¶

peagent_tool never bundles model weights. It resolves them, in order, from an explicit path argument, $PEAGENT_ASSET_ROOT, the local cache (~/.cache/peagent), and finally a Hugging Face download. describe_assets shows exactly what will be used, without loading anything.

In [2]:

from peagent_tool import describe_assets

pd.Series(describe_assets("maize"))

from peagent_tool import describe_assets pd.Series(describe_assets("maize"))

Out[2]:

species                                                           maize
asset_root                                                         None
cache_dir                                  /home/jy16611/.cache/peagent
hf_repo                                              YAOJ-bioin/peagent
hf_revision                                                        main
model_relative                                    models/maize/model.h5
model_path            /home/jy16611/.cache/peagent/models/maize/mode...
model_available                                                    True
metadata_relative                     metadata/maize_group_metadata.csv
metadata_path         /home/jy16611/.cache/peagent/metadata/maize_gr...
metadata_available                                                 True
dtype: object

2. Which cell types can I ask about?¶

Predictions are made for tissue_celltype groups, the exact labels used in the PEAgent atlas and browser. A group name is <tissue>_<celltype>, and only groups passing the minimum-cell threshold are kept. Always take the target name from this table rather than typing it by hand.

In [3]:

from peagent_tool import list_valid_celltypes

celltypes = list_valid_celltypes("maize")
print("valid tissue_celltype groups:", len(celltypes))
celltypes.head(10)

from peagent_tool import list_valid_celltypes celltypes = list_valid_celltypes("maize") print("valid tissue_celltype groups:", len(celltypes)) celltypes.head(10)

valid tissue_celltype groups: 94

Out[3]:

	group_idx	tissue_celltype	tissue	celltype	n_cells
0	0	leaf_bulliform	leaf	bulliform	1006
1	1	crownRoot_procambial_phloem_SE_precursors	crownRoot	procambial_phloem_SE_precursors	471
2	2	axillaryBud_axillary_L1_layer	axillaryBud	axillary_L1_layer	670
3	3	axillaryBud_ground_meristem	axillaryBud	ground_meristem	286
4	4	root_xylem_precursors	root	xylem_precursors	184
5	5	leaf_bundle_sheath	leaf	bundle_sheath	1515
6	6	crownRoot_cortex	crownRoot	cortex	3195
7	7	root_mixed	root	mixed	916
8	8	tassel_spikelet_meristem	tassel	spikelet_meristem	1174
9	9	axillaryBud_axillary_meristem	axillaryBud	axillary_meristem	515

In [4]:

celltypes.groupby("tissue").agg(n_celltypes=("celltype", "size"), n_cells=("n_cells", "sum"))

celltypes.groupby("tissue").agg(n_celltypes=("celltype", "size"), n_cells=("n_cells", "sum"))

Out[4]:

	n_celltypes	n_cells
tissue
axillaryBud	8	3088
crownRoot	23	15729
ear	17	5505
leaf	17	8029
root	12	3902
tassel	17	5001

3. Load the sequences¶

The model input window is 1,344 bp. Longer sequences are centre-cropped and shorter ones are centre-padded with N, so you can pass a sequence of any length and see what was done to it in the output columns.

The example FASTA holds two real maize ACR windows: a crownRoot-trichoblast cell-type-specific ACR and a constitutive ACR open in 93.5% of cell types.

In [5]:

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
CONSTITUTIVE = "chr1_5046150_5046651"   # constitutive ACR

{name: len(sequence) for name, sequence in records.items()}

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 CONSTITUTIVE = "chr1_5046150_5046651" # constitutive ACR {name: len(sequence) for name, sequence in records.items()}

Out[5]:

{'chr1_284657363_284657864': 1344, 'chr1_5046150_5046651': 1344}

4. Predict¶

load_predictor builds the model once and keeps it in memory; use it whenever you score more than one sequence. mode="exact" is the canonical PEAgent readout: the model predicts every single cell, then probabilities are averaged within each tissue_celltype group. This is exactly what the PEAgent web portal reports.

In [6]:

from peagent_tool import load_predictor

predictor = load_predictor("maize", mode="exact")
predictions = predictor.predict(records, top_k=10)
predictions[["sequence_id", "rank", "tissue_celltype", "tissue", "celltype", "n_cells", "prediction"]].head(10)

from peagent_tool import load_predictor predictor = load_predictor("maize", mode="exact") predictions = predictor.predict(records, top_k=10) predictions[["sequence_id", "rank", "tissue_celltype", "tissue", "celltype", "n_cells", "prediction"]].head(10)

Out[6]:

	sequence_id	rank	tissue_celltype	tissue	celltype	n_cells	prediction
0	chr1_284657363_284657864	1	root_protophloem	root	protophloem	212	0.844895
1	chr1_284657363_284657864	2	ear_provascular_inflorescence_strands	ear	provascular_inflorescence_strands	212	0.186956
2	chr1_284657363_284657864	3	crownRoot_trichoblast	crownRoot	trichoblast	370	0.174131
3	chr1_284657363_284657864	4	crownRoot_cortex_endodermis_initials	crownRoot	cortex_endodermis_initials	407	0.169712
4	chr1_284657363_284657864	5	root_lateral_root_epidermis	root	lateral_root_epidermis	617	0.160555
5	chr1_284657363_284657864	6	crownRoot_cortex	crownRoot	cortex	3195	0.155302
6	chr1_284657363_284657864	7	crownRoot_columella/epidermal_initials	crownRoot	columella/epidermal_initials	318	0.148011
7	chr1_284657363_284657864	8	crownRoot_endodermis	crownRoot	endodermis	235	0.143764
8	chr1_284657363_284657864	9	root_pericycle	root	pericycle	362	0.134397
9	chr1_284657363_284657864	10	crownRoot_phloem_sieve_element_precursors	crownRoot	phloem_sieve_element_precursors	427	0.130303

Nine of the ten top-ranked cell types for the specific ACR are root cell types, from the two root tissues the atlas samples (root and crownRoot). The constitutive ACR ranks differently — compare it directly:

In [7]:

predictions.query("sequence_id == @CONSTITUTIVE")[
    ["rank", "tissue_celltype", "prediction"]
].head(10)

predictions.query("sequence_id == @CONSTITUTIVE")[ ["rank", "tissue_celltype", "prediction"] ].head(10)

Out[7]:

	rank	tissue_celltype	prediction
10	1	root_protophloem	0.980594
11	2	ear_provascular_inflorescence_strands	0.881654
12	3	axillaryBud_axillary_L1_layer	0.663236
13	4	axillaryBud_ground_meristem	0.653322
14	5	crownRoot_cortex_endodermis_initials	0.652491
15	6	tassel_glume_primordia	0.627577
16	7	leaf_L1_leaf_primordia	0.607205
17	8	leaf_guard_mother_cell	0.603441
18	9	ear_floral_meristem	0.600049
19	10	crownRoot_cortex	0.577472

5. The full 94-cell-type profile¶

predict_wide returns one row per sequence and one column per tissue_celltype, which is the convenient form for scoring libraries of synthetic sequences.

In [8]:

wide = predictor.predict_wide(records)
print(wide.shape)
wide.iloc[:, :6]

wide = predictor.predict_wide(records) print(wide.shape) wide.iloc[:, :6]

(2, 95)

Out[8]:

	sequence_id	leaf_bulliform	crownRoot_procambial_phloem_SE_precursors	axillaryBud_axillary_L1_layer	axillaryBud_ground_meristem	root_xylem_precursors
0	chr1_284657363_284657864	0.021664	0.077106	0.056430	0.078449	0.128357
1	chr1_5046150_5046651	0.195004	0.503744	0.663236	0.653322	0.331178

In [9]:

summary = pd.DataFrame(
    {
        "cell-type-specific ACR": wide.set_index("sequence_id").loc[CTACR],
        "constitutive ACR": wide.set_index("sequence_id").loc[CONSTITUTIVE],
    }
)
summary.describe().T[["mean", "min", "50%", "max"]]

summary = pd.DataFrame( { "cell-type-specific ACR": wide.set_index("sequence_id").loc[CTACR], "constitutive ACR": wide.set_index("sequence_id").loc[CONSTITUTIVE], } ) summary.describe().T[["mean", "min", "50%", "max"]]

Out[9]:

	mean	min	50%	max
cell-type-specific ACR	0.048270	0.001101	0.018364	0.844895
constitutive ACR	0.256358	0.011452	0.177189	0.980594

Plotting the profile makes the difference in specificity obvious. Cell types are ordered by their score for the specific ACR.

In [10]:

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": 12,
})

top = summary.sort_values("cell-type-specific ACR", ascending=False).head(20)[::-1]
fig, ax = plt.subplots(figsize=(7.5, 6.5))
y = np.arange(len(top))
ax.barh(y - 0.2, top["cell-type-specific ACR"], height=0.4, color="#D7191C", label="cell-type-specific ACR")
ax.barh(y + 0.2, top["constitutive ACR"], height=0.4, color="#2C7BB6", label="constitutive ACR")
ax.set_yticks(y)
ax.set_yticklabels(top.index, fontsize=10)
ax.set_xlabel("Predicted accessibility")
ax.set_xlim(0, 1)
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(loc="lower right", frameon=False)
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": 12, }) top = summary.sort_values("cell-type-specific ACR", ascending=False).head(20)[::-1] fig, ax = plt.subplots(figsize=(7.5, 6.5)) y = np.arange(len(top)) ax.barh(y - 0.2, top["cell-type-specific ACR"], height=0.4, color="#D7191C", label="cell-type-specific ACR") ax.barh(y + 0.2, top["constitutive ACR"], height=0.4, color="#2C7BB6", label="constitutive ACR") ax.set_yticks(y) ax.set_yticklabels(top.index, fontsize=10) ax.set_xlabel("Predicted accessibility") ax.set_xlim(0, 1) 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(loc="lower right", frameon=False) fig.tight_layout() plt.show()

6. Scoring a synthetic sequence¶

Nothing above requires the sequence to exist in a genome. To design and screen synthetic regulatory sequences, build the string yourself and pass it in. Here a random background is used as a control and a Dof-family core motif is tiled into its centre.

In [11]:

rng = np.random.default_rng(2026)
background = "".join(rng.choice(list("ACGT"), size=1344, p=[0.29, 0.21, 0.21, 0.29]))

def implant(sequence, motif, centre=672, copies=1, spacing=12):
    sequence = list(sequence)
    start = centre - (copies * spacing) // 2
    for copy in range(copies):
        offset = start + copy * spacing
        sequence[offset : offset + len(motif)] = list(motif)
    return "".join(sequence)

designs = {
    "background": background,
    "background + 1x AAAGC": implant(background, "AAAGC", copies=1),
    "background + 4x AAAGC": implant(background, "AAAGC", copies=4),
}
design_scores = predictor.predict_wide(designs).set_index("sequence_id")
design_scores[["crownRoot_trichoblast", "root_trichoblast", "leaf_mesophyll", "ear_epidermal"]]

rng = np.random.default_rng(2026) background = "".join(rng.choice(list("ACGT"), size=1344, p=[0.29, 0.21, 0.21, 0.29])) def implant(sequence, motif, centre=672, copies=1, spacing=12): sequence = list(sequence) start = centre - (copies * spacing) // 2 for copy in range(copies): offset = start + copy * spacing sequence[offset : offset + len(motif)] = list(motif) return "".join(sequence) designs = { "background": background, "background + 1x AAAGC": implant(background, "AAAGC", copies=1), "background + 4x AAAGC": implant(background, "AAAGC", copies=4), } design_scores = predictor.predict_wide(designs).set_index("sequence_id") design_scores[["crownRoot_trichoblast", "root_trichoblast", "leaf_mesophyll", "ear_epidermal"]]

Out[11]:

	crownRoot_trichoblast	root_trichoblast	leaf_mesophyll	ear_epidermal
sequence_id
background	0.004513	0.005664	0.002133	0.002035
background + 1x AAAGC	0.004936	0.005605	0.001955	0.001987
background + 4x AAAGC	0.005720	0.006516	0.001879	0.001737

7. exact vs fast¶

mode="fast" averages the final dense-layer parameters within each cell group and applies the sigmoid to the grouped logit, instead of predicting all cells and then averaging probabilities. It is an approximation, so use exact for published numbers and fast for large screens.

In [12]:

import time

fast = load_predictor("maize", mode="fast")

start = time.time(); exact_wide = predictor.predict_wide(records); exact_seconds = time.time() - start
start = time.time(); fast_wide = fast.predict_wide(records); fast_seconds = time.time() - start

exact_values = exact_wide.set_index("sequence_id").to_numpy().ravel()
fast_values = fast_wide.set_index("sequence_id").to_numpy().ravel()
pd.Series({
    "exact seconds": round(exact_seconds, 2),
    "fast seconds": round(fast_seconds, 2),
    "pearson r": round(float(np.corrcoef(exact_values, fast_values)[0, 1]), 4),
    "max abs difference": round(float(np.abs(exact_values - fast_values).max()), 4),
})

import time fast = load_predictor("maize", mode="fast") start = time.time(); exact_wide = predictor.predict_wide(records); exact_seconds = time.time() - start start = time.time(); fast_wide = fast.predict_wide(records); fast_seconds = time.time() - start exact_values = exact_wide.set_index("sequence_id").to_numpy().ravel() fast_values = fast_wide.set_index("sequence_id").to_numpy().ravel() pd.Series({ "exact seconds": round(exact_seconds, 2), "fast seconds": round(fast_seconds, 2), "pearson r": round(float(np.corrcoef(exact_values, fast_values)[0, 1]), 4), "max abs difference": round(float(np.abs(exact_values - fast_values).max()), 4), })

Out[12]:

exact seconds         0.0200
fast seconds          0.4700
pearson r             0.9909
max abs difference    0.1163
dtype: float64

8. Sequence normalization¶

Any length is accepted. The returned table records the original length, the normalized length and what was done.

In [13]:

short = predictor.predict({"short_600bp": records[CTACR][372:972]}, top_k=3)
short[["sequence_id", "rank", "tissue_celltype", "prediction",
       "original_length", "normalized_length", "normalization_action"]]

short = predictor.predict({"short_600bp": records[CTACR][372:972]}, top_k=3) short[["sequence_id", "rank", "tissue_celltype", "prediction", "original_length", "normalized_length", "normalization_action"]]

Out[13]:

	sequence_id	rank	tissue_celltype	prediction	original_length	normalized_length	normalization_action
0	short_600bp	1	root_protophloem	0.826947	600	1344	center_padded_to_1344bp
1	short_600bp	2	crownRoot_trichoblast	0.161313	600	1344	center_padded_to_1344bp
2	short_600bp	3	crownRoot_cortex_endodermis_initials	0.159492	600	1344	center_padded_to_1344bp

9. The same thing from the command line¶

peagent-tool list-celltypes --species maize --out maize_celltypes.tsv

peagent-tool predict-celltypes \
  --species maize \
  --fasta maize_example_acrs.fa \
  --mode exact \
  --top-k 10 \
  --out predictions.tsv

Next¶

• Contribution maps — which bases in the sequence drive the prediction, globally and per cell type.
• Models and data — how to obtain the model weights and metadata.