Attribution and Motif Analysis

This tutorial covers per-base attribution, seqlet extraction, and TF-MoDISco motif discovery using a trained Cherimoya model.

For ref/alt and per-variant scoring see Variant Effect Prediction.

What are attributions?

Attributions quantify how much each base in the input sequence contributes to the model’s prediction. Cherimoya uses saturation mutagenesis — evaluating every single-nucleotide substitution and taking the predicted delta — to compute hypothetical importance scores. The output array has shape (n_examples, 4, window): one score per base per position.

These scores are commonly used to:

  • identify transcription factor binding sites,

  • discover de novo motifs (via TF-MoDISco), and

  • understand regulatory grammar.

Computing attributions (CLI)

cherimoya attribute -p attribute_params.json

Example JSON:

{
    "model": "my_model.torch",
    "sequences": "hg38.fa",
    "loci": "peaks.narrowPeak",
    "chroms": ["chr2", "chr4", "chr5"],
    "output": "counts",
    "batch_size": 512,
    "device": "cuda",
    "ohe_filename": "attributions.ohe.npz",
    "attr_filename": "attributions.attr.npz",
    "idx_filename": "attributions.idx.npy"
}

output controls what is being attributed to:

  • "counts" — attribute to total predicted log-counts (recommended for most analyses; uses cherimoya.LogCountWrapper).

  • "profile" — attribute to the predicted profile shape (uses cherimoya.ProfileWrapper).

The CLI automatically:

  1. Loads sequences from loci on chroms and filters out any example containing an N over the input window.

  2. Wraps the model with cherimoya.ControlWrapper (passing zero controls if the model has none) and then with the chosen output wrapper.

  3. Runs tangermeme.saturation_mutagenesis.saturation_mutagenesis over the central 400 bp of each input.

  4. Writes one-hot encoded inputs to ohe_filename, hypothetical importance scores to attr_filename, and a boolean mask (idx_filename) recording which loci survived the N-filter, so that downstream stages can re-align the attribution rows back to the original locus list.

The .npz outputs store the array under key arr_0 (numpy.savez_compressed default).

Computing attributions (Python)

import torch
from cherimoya import Cherimoya
from cherimoya import ControlWrapper
from cherimoya import LogCountWrapper
from cherimoya import ProfileWrapper
from tangermeme.saturation_mutagenesis import saturation_mutagenesis

model = Cherimoya.load("my_model.torch", device="cuda")

# ControlWrapper wraps the model so that .forward(X) returns just the
# profile/counts tuple, supplying zero controls if the model has none.
model = ControlWrapper(model)
wrapper = LogCountWrapper(model)   # use ProfileWrapper(...) to attribute to profile shape

# ISM over the central 400 bp of each input sequence.
mid = X.shape[-1] // 2
X_attr = saturation_mutagenesis(
    wrapper, X,
    batch_size=512,
    device="cuda",
    hypothetical=True,
    start=mid - 200, end=mid + 200,
)

This produces hypothetical importance scores. To get actual importance, multiply elementwise by the one-hot encoding and sum across the channel axis:

importance = (X_attr * X[:, :, mid - 200:mid + 200]).sum(dim=1)

Identifying seqlets

Seqlets are contiguous subsequences with high attribution scores that likely correspond to functional elements — binding motifs, in practice. They are extracted via TF-MoDISco-style recursive seqlet calling on the (attribution × one-hot) signal.

CLI:

cherimoya seqlets -p seqlet_params.json

Python:

from tangermeme.seqlet import recursive_seqlets

importance = (X_attr * X_ohe).sum(dim=1)

seqlets = recursive_seqlets(
    importance,
    threshold=0.01,
    min_seqlet_len=4,
    max_seqlet_len=25,
    additional_flanks=3,
)

The default seqlet parameters mirror the CLI defaults (see CLI Reference). After the recursive call, the CLI converts example-relative coordinates to genome coordinates using tangermeme.utils.example_to_fasta_coords and writes a BED file.

tomtom-lite annotation

When the pipeline JSON provides a motifs MEME file, the pipeline subcommand additionally invokes ttl (tomtom-lite) on the seqlet BED to annotate each seqlet with its closest match against the motif database. This is what produces {name}.seqlets_annotated.bed and {name}.motif_seqlet_count.tsv.

If you want to run this independently, the equivalent shell call is:

ttl -f hg38.fa -b seqlets.bed \
    -t JASPAR_2024.meme \
    -s 100 -m 1000 -a 100 -c 250 -j -1 > seqlets_annotated.bed

TF-MoDISco motif discovery

TF-MoDISco clusters seqlets into motif patterns. The pipeline runs this automatically; run it manually like so:

modisco motifs \
    -s attributions.ohe.npz \
    -a attributions.attr.npz \
    -n 100000 \
    -o modisco_results.h5

modisco report \
    -i modisco_results.h5 \
    -o modisco_report/ \
    -s ./

The pipeline’s TF-MoDISco step uses 100,000 seqlets by default (modisco_motifs_parameters.n_seqlets).

Motif marginalization

To quantify the causal effect of an inserted motif on the predicted profile and counts:

cherimoya marginalize -p marginalize_params.json

The output directory contains per-motif plots and a summary report showing how predictions change when each motif is inserted at the center of negative (non-peak) backgrounds.

Note

Marginalization requires a motif database in MEME format. JASPAR provides such files for many species; the latest as of writing is JASPAR 2024.