15. Cross-species prediction¶

Bolero reads DNA (plus a cell-state embedding) — nothing about it is human- or mouse-specific. So the same trained atlas checkpoint runs on any species' genome, with no retraining. That is the whole idea behind Figure 6: predict cell-type-specific accessibility across mammals by pointing Bolero at each species' own genome, then compare the predictions through orthology.

This page demos the core function — running the trained atlas model on a non-human genome (macaque, rheMac10) in DNA-only mode. It is concise and demo-only: the cross-species comparison (peak liftover, conservation of accessibility / AP-1 activity) is a downstream analysis and out of scope here — this is the per-species prediction call that feeds it.

This uses the bolerodata model zoo and data lake (the paper route, like the variant-effect and score-model pages). The atlas model was trained across six mammals, so macaque is one of its datasets; the same call works for any species' genome — including one the model never saw (the mouse-model-predicts-human transfer in Fig 6b is this same call with a held-out genome).

Setup¶

In [1]:

from bolero import init
from bolerodata import MODELS

init(num_cpus=24)

from bolero import init from bolerodata import MODELS init(num_cpus=24)

2026-07-17 13:19:07,048	INFO worker.py:1781 -- Started a local Ray instance.

Pick the model and the target species¶

MODELS["Borzoi10M"] is the multi-dataset atlas ATAC model (the same one behind the accessibility page). Because it is multi-dataset you must tell it whose cell states you are conditioning on with model.dataset — here the macaque cortical cell states, since we predict on the macaque genome. Everything about the target species comes off its Dataset record: its genome assembly and its peak set.

In [2]:

model = MODELS["Borzoi10M"]
model.dataset = (
    "MultiSpeciesMacaque"  # multi-dataset model: select the conditioning cell states
)
target = model.dataset  # the Dataset object for the species we predict on

print("target species :", target.species)
print("target genome  :", target.genome)  # rheMac10 — a non-human assembly
print("target peaks   :", target.peak_bed_path)

model = MODELS["Borzoi10M"] model.dataset = ( "MultiSpeciesMacaque" # multi-dataset model: select the conditioning cell states ) target = model.dataset # the Dataset object for the species we predict on print("target species :", target.species) print("target genome :", target.genome) # rheMac10 — a non-human assembly print("target peaks :", target.peak_bed_path)

target species : Macaque
target genome  : rheMac10
target peaks   : /large_storage/zhoulab/hanliu/wmb/MultiSpeciesM1/MultiSpeciesM1.macaque.peaks.bed

Cross-species = swap the genome + run DNA-only¶

A within-species run and a cross-species run are the same task; only two things change.

1. The sequence source. Point the predictor at the target species' genome and peak_bed_path. Windows and peaks are then in that species' own coordinates (rheMac10), not lifted from human.
2. DNA-only inference (nosignal=True). A non-reference species has no reference-coverage track aligned to the model, so the prediction is driven purely by DNA + the cell-state embedding. (use_ref_bw=False drops the reference bigWig; embedding_only_mode=True skips ground-truth loading and metrics — there is no matched measured track to score against here.)

sample_n_pseudobulks=8 keeps just a handful of macaque cell states so the demo is quick.

In [3]:

predictor = model.create_predictor(
    subset_name=None,
    sample_n_pseudobulks=8,  # a handful of macaque cell states for the demo
    nosignal=True,  # DNA-only: the cross-species inference mode
    use_ref_bw=False,  # no reference bigWig (paired with nosignal)
    embedding_only_mode=True,  # drive from embeddings; skip truth/metrics
    genome=target.genome,  # <- swap the sequence source to macaque (rheMac10)
    peak_bed_path=target.peak_bed_path,
)

predictor = model.create_predictor( subset_name=None, sample_n_pseudobulks=8, # a handful of macaque cell states for the demo nosignal=True, # DNA-only: the cross-species inference mode use_ref_bw=False, # no reference bigWig (paired with nosignal) embedding_only_mode=True, # drive from embeddings; skip truth/metrics genome=target.genome, # <- swap the sequence source to macaque (rheMac10) peak_bed_path=target.peak_bed_path, )

Tile the target genome into windows¶

The model predicts on 524,288 bp windows. make_windows builds them straight from the macaque genome — so these are rheMac10 coordinates. We keep a handful for the demo; drop the slice to predict the whole genome (~5.7k windows), which is what Fig 6 does per species.

In [4]:

regions = predictor.borzoi_regions.genome.make_windows(
    524288, 500000, force_length=True, as_df=True
)
regions["Name"] = (
    regions["Chromosome"].astype(str)
    + ":"
    + regions["Start"].astype(str)
    + "-"
    + regions["End"].astype(str)
)

regions_demo = regions.iloc[
    :8
].copy()  # demo: a few macaque windows (drop to predict all)
print("macaque genome windows:", regions.shape[0], "-> predicting", len(regions_demo))
regions_demo.head()

regions = predictor.borzoi_regions.genome.make_windows( 524288, 500000, force_length=True, as_df=True ) regions["Name"] = ( regions["Chromosome"].astype(str) + ":" + regions["Start"].astype(str) + "-" + regions["End"].astype(str) ) regions_demo = regions.iloc[ :8 ].copy() # demo: a few macaque windows (drop to predict all) print("macaque genome windows:", regions.shape[0], "-> predicting", len(regions_demo)) regions_demo.head()

macaque genome windows: 5697 -> predicting 8

Out[4]:

	Chromosome	Start	End	Name
0	chr1	0	524288	chr1:0-524288
1	chr1	500000	1024288	chr1:500000-1024288
2	chr1	1000000	1524288	chr1:1000000-1524288
3	chr1	1500000	2024288	chr1:1500000-2024288
4	chr1	2000000	2524288	chr1:2000000-2524288

Run the prediction task¶

Same prediction_task as every other prediction page — we just hand it the macaque windows. It streams each window × macaque cell state through the model and writes the per-peak predictions to output_dir.

In [5]:

output_dir = "cross_species_macaque"
predictor.prediction_task(
    output_dir,
    regions=regions_demo,  # macaque (rheMac10) windows; drop for whole-genome
    mode="prediction",
)
print("wrote predictions to", output_dir)

output_dir = "cross_species_macaque" predictor.prediction_task( output_dir, regions=regions_demo, # macaque (rheMac10) windows; drop for whole-genome mode="prediction", ) print("wrote predictions to", output_dir)

No parquet dataset found, skipping region filtering.
Saving batches to /large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/docs/tutorials/cross_species/cross_species_macaque/batch
Saving stats to /large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/docs/tutorials/cross_species/cross_species_macaque/summary_stats.joblib.gz
Using temporary directory /tmp/tmpwezmqrpm
0 regions has finished in /large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/docs/tutorials/cross_species/cross_species_macaque/batch
8 regions to compute
Data loader batch size 16

Loading base model weights from: /home/hanliu/data/wmb/model/Borzoi/flashzoi_checkpoints/fold_0.renamed.pt

Model does not have gene_count_output_head, skip
Loading checkpoint from /large_storage/zhoulab/hanliu/250907-MultiDatasetModel/Bolero10M/model/251103-Bolero10M_atac_fold0-final.lora.best_checkpoint.pt

==========
Dataloader Batch Schema:
- MultiSpeciesMacaque.MetaCell:cov_scale: <class 'torch.Tensor'> torch.Size([16]) torch.float64 cuda:0
- __conditionemb__: <class 'torch.Tensor'> torch.Size([16, 1, 4]) torch.float32 cuda:0
- __conditionemb__:cond0: <class 'torch.Tensor'> torch.Size([8, 1, 4]) torch.float32 cuda:0
- __conditionemb__:cond1: <class 'torch.Tensor'> torch.Size([8, 1, 4]) torch.float32 cuda:0
- __conditionemb__:delta: <class 'torch.Tensor'> torch.Size([8, 1, 4]) torch.float32 cuda:0
- __dna__: <class 'torch.Tensor'> torch.Size([8, 4, 524288]) torch.uint8 cuda:0
- __embedding__: <class 'torch.Tensor'> torch.Size([16, 60]) torch.float32 cuda:0
- __embedding__:cond0: <class 'torch.Tensor'> torch.Size([8, 60]) torch.float32 cuda:0
- __embedding__:cond1: <class 'torch.Tensor'> torch.Size([8, 60]) torch.float32 cuda:0
- __embedding__:delta: <class 'torch.Tensor'> torch.Size([8, 60]) torch.float32 cuda:0
- __ypred__:cond1: <class 'torch.Tensor'> torch.Size([8, 8, 16352]) torch.float32 cuda:0
- __ypred__:peak:cond1: <class 'torch.Tensor'> torch.Size([8, 770]) torch.float32 cuda:0
- condition_pairs: <class 'numpy.ndarray'> (8, 2) <U8
- peak: <class 'pandas.core.frame.DataFrame'> (770, 7)
- pseudobulk_ids: <class 'numpy.ndarray'> (16,) <U24
- region: <class 'numpy.ndarray'> (8,) <U20
- region_name: <class 'numpy.ndarray'> (8,) <U20
==========

Total time: 5.39s
Inference time: 4.93s or 4.929s per batch
Callback time: 0.27s or 0.270s per batch
(total 1 batches)

==========
Saved Batch Schema:
- __embedding__: <class 'numpy.ndarray'> (16, 60) float32
- __ypred__:peak:cond1: <class 'numpy.ndarray'> (8, 770) float32
- condition_pairs: <class 'numpy.ndarray'> (8, 2) <U8
- peak: <class 'pandas.core.frame.DataFrame'> (770, 7)
- pseudobulk_ids: <class 'numpy.ndarray'> (16,) <U24
- region: <class 'numpy.ndarray'> (8,) <U20
- region_name: <class 'numpy.ndarray'> (8,) <U20
==========

==========
Final Stats Batch Schema:
- condition_pairs: <class 'numpy.ndarray'> (8, 2) <U8
- peak: <class 'pandas.core.frame.DataFrame'> (770, 7)
- pseudobulk_ids: <class 'numpy.ndarray'> (16,) <U24
- region: <class 'numpy.ndarray'> (8,) <U20
==========

Removed temporary files in /tmp/tmpwezmqrpm
1 batches in /large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/docs/tutorials/cross_species/cross_species_macaque
True peak data not found in batch, skip true peak data saving.
wrote predictions to cross_species_macaque

What you get¶

prediction_task writes pred_peak_data.feather — predicted per-peak accessibility, a peak × cell-state matrix (macaque peaks as rows, the sampled cell states as columns), all in rheMac10 coordinates. There is no measured track or accuracy metric here (DNA-only, embedding_only_mode); the output is purely the model's prediction on the target species' sequence.

In [6]:

import os

import pandas as pd

print("output files:", sorted(os.listdir(output_dir)))

pred = pd.read_feather(f"{output_dir}/pred_peak_data.feather")
print("\npred_peak_data.feather:", pred.shape, "(macaque peaks x cell states)")
pred.head()

import os import pandas as pd print("output files:", sorted(os.listdir(output_dir))) pred = pd.read_feather(f"{output_dir}/pred_peak_data.feather") print("\npred_peak_data.feather:", pred.shape, "(macaque peaks x cell states)") pred.head()

output files: ['batch', 'config.joblib.gz', 'pred_peak_data.feather', 'summary_stats.joblib.gz']

pred_peak_data.feather: (751, 8) (macaque peaks x cell states)

Out[6]:

	group0-pseudobulk59	group0-pseudobulk65	group0-pseudobulk78	group1-pseudobulk101	group1-pseudobulk165	group1-pseudobulk193	group1-pseudobulk69	group1-pseudobulk70
chr1:101596-102096	26.491409	30.271399	27.159431	11.286598	12.375208	12.141721	10.887924	13.884491
chr1:102320-102820	153.394836	191.166519	138.662033	75.595825	70.910721	68.313652	65.677383	100.567032
chr1:102975-103475	16.162777	15.712399	17.185543	13.270825	13.888287	12.820403	11.658803	13.981735
chr1:104684-105184	2.323617	2.381340	3.577239	12.104186	7.872944	6.501991	9.217743	3.233615
chr1:105703-106203	2.725699	2.976972	3.403375	10.311110	8.685028	8.263292	8.611458	6.409596

From one species to a phylogeny¶

That is the entire cross-species function: point the model at a species' genome + peaks and run DNA-only. Fig 6 repeats this call per species — human (hg38), macaque (rheMac10), marmoset (calJac4), mouse (mm10), and outward to 85 Zoonomia mammals — then aligns the per-species predictions through peak liftover / orthology to ask which cell-type accessibility (and which AP-1 activity) is conserved. That comparison is the downstream analysis; this page is the prediction step it stands on.

For the per-species AP-1 activity panels of Fig 6, swap Borzoi10M for the AP-1 score model Bolero10M+AP1+260104 (see the score-model prediction page) — the cross-species call is identical, each cell state just also carries its AP-1 score.