Tutorial: Single-Cell-Resolution Sequence Models with Scooby using cellink

This tutorial demonstrates how to use scooby (Hingerl et al. 2025) through the cellink package to predict single-cell-resolution scRNA-seq coverage directly from DNA sequence, and to score the effect of individual variants on that coverage.

Scooby fine-tunes the pretrained multi-omics profile predictor Borzoi with a small cell-specific decoder (via LoRA), conditioned on a precomputed per-cell embedding. Unlike bulk sequence-to-expression models, it predicts coverage for a specific cell (or a pseudobulk of cells) rather than a fixed panel of bulk tracks.

Everything in this notebook runs in a single environment via pip install cellink[scooby], plus cellink[embpy] for the one variant-scoring cell. The one exception, an alternative, more faithful embedding-building method (build_scooby_embedding_scpoli), needs a separate environment, so it’s discussed, rather than run here.

pip install cellink[scooby]   # torch, accelerate, enformer-pytorch, borzoi-pytorch, scooby itself
pip install cellink[embpy]    # only needed for the variant-scoring cell near the end

import numpy as np
import pandas as pd

from cellink.resources import get_dummy_onek1k
from cellink.tl.external import (
    build_scooby_embedding,
    configure_scooby_runner,
    load_scooby_checkpoint,
    predict_scooby_profile,
    KNOWN_SCOOBY_CHECKPOINTS,
)

print("Released checkpoints this integration knows about:")
for name, info in KNOWN_SCOOBY_CHECKPOINTS.items():
    print(f"  {name}: {info}")

[2026-07-25 11:13:23,699] WARNING:cellink.resources._datasets_utils: liftover unavailable (No module named 'liftover'); hg19<->hg38 coordinate translation will be disabled.
Released checkpoints this integration knows about:
  johahi/neurips-scooby: {'cell_emb_dim': 14, 'n_tracks': 3, 'modality': 'multiome', 'use_transform_borzoi_emb': True}
  lauradmartens/onek1k-scooby: {'cell_emb_dim': 10, 'n_tracks': 2, 'modality': 'rna', 'use_transform_borzoi_emb': True}
  lauradmartens/epicardioids-scooby: {'cell_emb_dim': 50, 'n_tracks': 3, 'modality': 'multiome', 'use_transform_borzoi_emb': True}

Loading a Released Checkpoint and Predicting Coverage

load_scooby_checkpoint resolves cell_emb_dim/n_tracks automatically for known released checkpoints (see KNOWN_SCOOBY_CHECKPOINTS above), for your own fine-tunes, pass them explicitly (they aren’t recoverable from the checkpoint’s own config.json, which only carries Borzoi’s base hyperparameters).

We fetch a real 524,288bp genomic window around GAPDH (a near-universally expressed housekeeping gene, GRCh38 chr12) via embpy’s SequenceProvider, which falls back to the Ensembl REST API when no local reference FASTA is given.

configure_scooby_runner(device="auto")
model = load_scooby_checkpoint("lauradmartens/onek1k-scooby")
print(f"Loaded {model.__class__.__name__}, cell_emb_dim={model.cell_emb_dim}, n_tracks={model.n_tracks}")

from embpy.tl.genomics import SequenceProvider

provider = SequenceProvider()  # no fasta_file given -> Ensembl REST fallback
window, offset = provider.get_window("chr12", 6_536_000, context=524_288)
print(f"Fetched a real {len(window)}bp window around GAPDH; SNP offset within window: {offset}")

base_to_idx = {"A": 0, "C": 1, "G": 2, "T": 3}
one_hot_seq = np.zeros((len(window), 4), dtype=np.float32)
for i, b in enumerate(window.upper()):
    if b in base_to_idx:
        one_hot_seq[i, base_to_idx[b]] = 1.0

# 10 cell embeddings (cell_emb_dim=10 for this checkpoint) 
# using random vectors here purely to keep this cell self-contained.
rng = np.random.default_rng(0)
cell_embs = rng.normal(scale=0.3, size=(10, model.cell_emb_dim)).astype(np.float32)

profile = predict_scooby_profile(model, one_hot_seq, cell_embs, aggregate="pseudobulk")
print(f"Predicted pseudobulk profile: {profile.shape} (n_tracks x n_bins, 32bp bins)")
print(f"Max predicted coverage: {profile.max():.2f} (RNA:+ and RNA:- strands)")

Loaded Scooby, cell_emb_dim=10, n_tracks=2
[2026-07-25 11:14:33,516] INFO:rdkit: Enabling RDKit 2026.03.4 jupyter extensions
[2026-07-25 11:14:34,638] INFO:root: SequenceProvider: REST fetch https://rest.ensembl.org/sequence/region/human/12:6273856..6798143:1?content-type=text/plain ...
Fetched a real 524288bp window around GAPDH; SNP offset within window: 262145
Predicted pseudobulk profile: (6144, 2) (n_tracks x n_bins, 32bp bins)
Max predicted coverage: 12.99 (RNA:+ and RNA:- strands)

Scoring a Variant’s Effect

For variant-effect scoring (diffing predicted coverage between the reference and alternate allele) we utilize score_variant_effects_scooby, which utilizes the logic implemented in the embpy package (pip install cellink[embpy]).

from cellink.tl.external import score_variant_effects_scooby
from embpy.tl.genomics import SNPContext

snp_pos = offset
ref_allele = window[snp_pos - 1].upper() 
alt_allele = "G" if ref_allele != "G" else "T"
snp = SNPContext(
    chrom="chr12",
    position=snp_pos, 
    ref_allele=ref_allele,
    alt_alleles=[alt_allele],
    context_window=len(window),
    strand="+",
    variant_id="demo_variant",
)

vep = score_variant_effects_scooby(
    snp, window, "lauradmartens/onek1k-scooby", cell_embs, aggregate="pseudobulk",
)
print(f"Effect scores shape: {vep.effect_scores[0].shape}")
print(f"Max |log2FC| ({ref_allele}->{alt_allele} at window offset {snp_pos}): {np.max(np.abs(vep.effect_scores[0])):.4f}")

[2026-07-25 11:16:03,842] INFO:root: Loading Scooby 'lauradmartens/onek1k-scooby' …
[2026-07-25 11:16:04,163] INFO:root: Scooby 'lauradmartens/onek1k-scooby' loaded on cpu (cell_emb_dim=10, n_tracks=2).
Effect scores shape: (2,)
Max |log2FC| (A->G at window offset 262145): 0.0035

Load Data

We use a dummy OneK1K dataset here (~100 donors, real single-cell expression, real cell-type labels). It is small enough to build an embedding interactively.

dd = get_dummy_onek1k()
print(f"Cells: {dd.C.shape[0]}, Genes: {dd.C.shape[1]}")
print("Available obs columns:", dd.C.obs.columns.tolist())

[2026-07-25 11:13:42,510] INFO:root: /nfs/users/nfs_l/la17/cellink_data/dummy_onek1k/dummy_onek1k.dd.h5 already exists
[2026-07-25 11:13:42,511] INFO:root: Veryifying checksum
[2026-07-25 11:13:51,571] INFO:root: Loaded dummy OneK1K dataset: (100, 146939, 125366, 34073)
Cells: 125366, Genes: 34073
Available obs columns: ['orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'donor_id', 'pool_number', 'predicted.celltype.l2', 'predicted.celltype.l2.score', 'age', 'organism_ontology_term_id', 'tissue_ontology_term_id', 'assay_ontology_term_id', 'disease_ontology_term_id', 'cell_type_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'development_stage_ontology_term_id', 'sex_ontology_term_id', 'is_primary_data', 'suspension_type', 'tissue_type', 'cell_type', 'assay', 'disease', 'organism', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid']

Building a Per-Cell Embedding

To fine-tune your own checkpoint (see the next section) you first need a per-cell embedding to condition on. Any (n_cells, cell_emb_dim) array works as input to scooby. build_scooby_embedding is the quickest way to get one. It expects raw counts in adata.X.

This embedding is appropriate for training your own checkpoint from scratch, but its scale depends entirely on your own data, so it will not match what a released checkpoint (like the one used above) was actually trained against. cellink also provides build_scooby_embedding_scpoli, a more faithful recipe based on scPoli that matches the real embedding released checkpoints were trained with. It needs a separate Python environment, so it is discussed in theory below.

embedding = build_scooby_embedding(dd.C.copy(), n_comps=16)
emb = np.stack(embedding["embedding"].to_numpy())
print(f"Embedding: {emb.shape}, mean={emb.mean():.3f}, std={emb.std():.3f}")
embedding.head()

Embedding: (125366, 16), mean=-0.000, std=1.887

	embedding
barcode
ATGAGGGAGTACGCCC-15	[6.483895, 0.019763561, 0.8057307, 2.2672951, ...
AGCAGCCGTTACTGAC-15	[7.861376, 0.86356413, -0.49554238, -2.2540953...
AGCATACAGATGTAAC-15	[2.6045034, 0.84194547, 2.403064, 3.0084584, 1...
AGCCTAACAAACGCGA-15	[2.4460804, 0.07007817, 0.8245969, 3.5522282, ...
CGTGTCTGTGATGTGG-15	[-3.0054328, -2.6365871, -1.1180134, -0.754107...

build_scooby_embedding_scpoli matches the real recipe (scPoli, negative-binomial reconstruction loss, seurat_v3-flavor HVG selection, conditioned on batch/sample and cell type). It genuinely needs a separate Python environment: scarches (0.6.x, the version this integration was built against) pins conflict with the newer scvi-tools/torch versions scooby itself needs. Install cellink[scpoli] in its own dedicated environment, not alongside cellink[scooby].

pip install cellink[scpoli]   # scarches + scvi-tools, in a separate env from cellink[scooby]

from cellink.resources import get_dummy_onek1k
from cellink.tl.external import build_scooby_embedding_scpoli
import numpy as np

dd = get_dummy_onek1k()
dd.aggregate(obs=["donor_id"], func="first", add_to_obs=True)

embedding_scpoli = build_scooby_embedding_scpoli(
    dd.C.copy(),
    condition_key="donor_id",
    cell_type_key="predicted.celltype.l2",
    latent_dim=16,
    n_epochs=5,
    pretraining_epochs=5,
    n_top_genes=500,
    checkpoint_dir="./scpoli_tutorial_checkpoint",
)
emb_scpoli = np.stack(embedding_scpoli["embedding"].to_numpy())
print(f"scPoli embedding: {emb_scpoli.shape}, mean={emb_scpoli.mean():.3f}, std={emb_scpoli.std():.3f}")

Fine-Tuning Your Own Checkpoint

train_scooby (RNA-only) and train_scooby_multiome (RNA+ATAC) wrap the same LoRA fine-tuning recipe Scooby’s own reference training scripts use. This cell is deliberately not executed here: a real training run needs real fragment-coverage h5ads, a real genome FASTA, and hours of GPU time even for a small pilot-scale run. The call itself looks like:

from cellink.tl.external import train_scooby

train_scooby(
    rna_plus_path="your_data_plus.h5ad",
    rna_minus_path="your_data_minus.h5ad",
    embedding_path="embedding.pq",       # e.g. written by build_scooby_embedding above
    output_dir="checkpoints",
    run_name="my_fine_tune",
    sequences_path="sequences_human.bed",
    genome_path="genome.fa",
    cell_emb_dim=16,                     # must match your embedding's dimensionality
)

This trains via accelerate, which checkpoints with accelerator.save_state(), a model.safetensors with no config.json, meant for resuming training, not for Scooby.from_pretrained(). To get a checkpoint you can actually load for inference (as load_scooby_checkpoint does above), convert it first:

from cellink.tl.external import convert_scooby_lora_checkpoint

convert_scooby_lora_checkpoint(
    checkpoint_dir="checkpoints/my_fine_tune/final",
    output_dir="checkpoints/my_fine_tune/final_pretrained",
    cell_emb_dim=16,
    n_tracks=2,        # 2 for RNA-only (plus/minus strand), 3 for multiome (+ ATAC)
)

The result can be passed directly to load_scooby_checkpoint/predict_scooby_profile/score_variant_effects_scooby exactly like the released checkpoint used above.