Motifs and deep learning¶
Sequence models such as chromBPNet learn to predict a genomic signal, here chromatin accessibility, from the DNA sequence alone. What a model has learned can be read back by asking how much each base of an input contributed to its prediction, with attribution methods such as DeepLIFT, and TF-MoDISco then gathers the short stretches of high contribution that recur across thousands of regions into patterns, much as a motif discovery tool would. This guide works through the published outputs of one such model, from its TF-MoDISco patterns down to the contribution scores of a single enhancer, and shows how motifs and sequences pass between the package and a model.
The model is the chromBPNet model of ATAC-seq in the erythroleukemia cell
line K562 that ENCODE publishes as annotation
ENCSR893SUD.
Nothing here trains or runs a model, since everything we need was written
out when the model was interpreted, so no deep-learning library is needed,
though the arrays we make can be passed to PyTorch as they are. Reading
TF-MoDISco results needs the h5py extra and reading bigWig files the
bigwig extra.
The guide needs version 0.2.0 or later, which until its release means
installing the development version.
import json
import pathlib
import tarfile
import urllib.request
import matplotlib.pyplot as plt
import numpy as np
import polars as pl
import universalmotif as um
pl.Config.set_tbl_rows(25)
pl.Config.set_fmt_str_lengths(40);
Downloading the data¶
We need two files from ENCODE, the TF-MoDISco results of the model and the list of regions it was interpreted on, along with the JASPAR 2024 collection of vertebrate motifs. As in the other guides, a small helper downloads each file once into a cache directory in your home folder. The TF-MoDISco results come as a 311 MB archive holding the patterns found for both of the model's outputs. We only use those of the counts output, which predicts how accessible each region is, so only that file is unpacked:
CACHE = pathlib.Path.home() / ".cache" / "universalmotif-docs"
def fetch(url, name=None):
CACHE.mkdir(parents=True, exist_ok=True)
path = CACHE / (name or url.rsplit("/", 1)[-1])
if not path.exists():
part = path.with_name(path.name + ".part")
urllib.request.urlretrieve(url, part)
part.rename(path)
return path
ENCODE = "https://www.encodeproject.org/files"
modisco_tar = fetch(f"{ENCODE}/ENCFF009WUL/@@download/ENCFF009WUL.tar.gz")
regions_path = fetch(f"{ENCODE}/ENCFF834JBL/@@download/ENCFF834JBL.bed.gz")
jaspar_path = fetch("https://jaspar.elixir.no/download/data/2024/CORE/JASPAR2024_CORE_vertebrates_non-redundant_pfms_meme.txt")
modisco_path = CACHE / "counts" / "tfmodisco.raw_output.counts.ENCSR868FGK.hd5"
if not modisco_path.exists():
with tarfile.open(modisco_tar) as tar:
tar.extract(f"counts/{modisco_path.name}", CACHE, filter="data")
The contribution score of every base is in a third file, a bigWig track of 2.1 GB, which we don't download at all. A bigWig file is indexed, so the parts of the track we want can be read straight from ENCODE's server. The genome sequence of those parts comes from the UCSC Genome Browser's API, again without downloading the whole genome. Both count positions from 0 and leave the end out, as the package does:
scores_url = f"{ENCODE}/ENCFF407GCO/@@download/ENCFF407GCO.bigWig"
def hg38(chrom, start, end):
url = f"https://api.genome.ucsc.edu/getData/sequence?genome=hg38;chrom={chrom};start={start};end={end}"
with urllib.request.urlopen(url) as response:
return json.load(response)["dna"].upper()
TF-MoDISco patterns¶
TF-MoDISco writes its results to an HDF5 file holding, for each pattern,
the seqlets it was built from and their average contribution scores, the
pattern's contribution weight matrix or CWM.
read_modisco() reads every pattern as a
motif of type CWM, trimmed to the positions whose contribution is at least
30% of the highest:
cwms = um.read_modisco(modisco_path)
print(len(cwms))
cwms[0].name
60
'pos_patterns.ENCSR868FGK_ATAC_K562_ENCSR467RSV_counts_pattern_0'
The patterns come in two groups. Positive patterns are made of bases that raised the model's prediction, and negative patterns of bases that lowered it. Their names also carry the accessions of the experiment and the model, which we drop, keeping the group and the number:
def short_name(motif):
group, _, name = motif.name.partition(".")
return motif.replace(name=f"{group[:3]}_{name.rsplit('_', 1)[1]}")
cwms = [short_name(m) for m in cwms]
pos = [m for m in cwms if m.name.startswith("pos")]
neg = [m for m in cwms if m.name.startswith("neg")]
print(len(pos), len(neg))
pl.DataFrame({
"pattern": [m.name for m in cwms],
"seqlets": [int(m.nsites) for m in cwms],
"width": [m.width for m in cwms],
"consensus": [m.consensus for m in cwms],
}).head(8)
44 16
| pattern | seqlets | width | consensus |
|---|---|---|---|
| str | i64 | i64 | str |
| "pos_0" | 36546 | 5 | "GATAA" |
| "pos_1" | 33538 | 14 | "CCACYAGGGGGCGS" |
| "pos_2" | 27345 | 9 | "GCCCCGCCC" |
| "pos_3" | 13390 | 10 | "GCTGAGTCAY" |
| "pos_4" | 9410 | 7 | "CTTCCGG" |
| "pos_5" | 6024 | 7 | "YGATTGG" |
| "pos_6" | 4957 | 11 | "GCGCATGCGCR" |
| "pos_7" | 4175 | 8 | "TTGCATCA" |
Within each group the patterns are ordered by their number of seqlets, so
the most common pattern the model relies on is the short GATAA, found
over 36,000 times. A CWM gives the average contribution of each letter at
each position, so its logo shows how much each base mattered to the model
rather than how often it was seen, and the letters of negative patterns
hang below zero. Here are the first five positive patterns and the first
two negative ones, drawn as they are rather than aligned to one another:
fig = um.view_motifs(pos[:5] + neg[:2], use_type="CWM", rc=False, min_overlap=100,
names_pos="right", show_positions=False, figsize=(7, 8))
All seven are drawn on the same scale, which shows how much weaker the negative patterns are than the strongest positive ones.
Naming the patterns¶
TF-MoDISco finds patterns without knowing which proteins bind them, so the
next step is to match them to known motifs. Comparisons work on
probabilities, and with type="PPM", read_modisco() gives each pattern as
the letter frequencies of its seqlets instead. We also trim less here,
keeping positions with at least 20% of the highest contribution, for a
reason the first pattern shows:
ppms = [short_name(m) for m in um.read_modisco(modisco_path, type="PPM", trim_threshold=0.2)]
print(pos[0].consensus, ppms[0].consensus)
GATAA KNNNNNNNNWGATAAG
At the default cutoff the first pattern is just the five letters GATAA,
which is too short to stand out against a database of nearly 900 motifs,
whereas the wider version keeps enough of its flanks to be matched with
confidence. We compare every pattern to JASPAR with
compare_motifs(), as in the
ChIP-seq workflow, and keep the best match of
each:
jaspar = um.read_meme(jaspar_path)
tf_names = {m.name: m.altname for m in jaspar}
cmp = um.compare_motifs(ppms + jaspar, compare_to=list(range(len(ppms))), qvalue=0.05, nthreads=2)
best = (
cmp.filter(pl.col("target_i") >= len(ppms))
.with_columns(tf=pl.col("target").replace_strict(tf_names))
.sort("qvalue")
.group_by("subject", maintain_order=True)
.first()
)
summary = pl.DataFrame({
"pattern": [m.name for m in ppms],
"seqlets": [int(m.nsites) for m in ppms],
"consensus": [m.consensus for m in ppms],
}).join(best.select(pattern="subject", tf="tf", qvalue="qvalue"), on="pattern", how="left")
summary.head(20)
| pattern | seqlets | consensus | tf | qvalue |
|---|---|---|---|---|
| str | i64 | str | str | f64 |
| "pos_0" | 36546 | "KNNNNNNNNWGATAAG" | "GATA2" | 0.003661 |
| "pos_1" | 33538 | "CCASYAGRKGGCRSY" | "CTCF" | 1.7152e-14 |
| "pos_2" | 27345 | "NGCCCCGCCCCC" | "SP4" | 0.000058 |
| "pos_3" | 13390 | "GMTGAGTCAY" | "MAF::NFE2" | 0.000052 |
| "pos_4" | 9410 | "NRCTTCCGG" | "ELK4" | 0.000026 |
| "pos_5" | 6024 | "NNNYNCYSATTGGCYN" | "NFYA" | 0.010388 |
| "pos_6" | 4957 | "NGCGCNTGCGCN" | "Nrf1" | 0.000028 |
| "pos_7" | 4175 | "RTTGCATCAK" | "ATF4" | 0.000012 |
| "pos_8" | 3682 | "NCACTTCCTSTTT" | "SPIB" | 0.000024 |
| "pos_9" | 2785 | "WWNTGCWNWGTCAY" | "MAF::NFE2" | 0.001445 |
| "pos_10" | 2761 | "NNNNNNNNNRTGACGYMA" | "FOSL1::JUND" | 0.000902 |
| "pos_11" | 2623 | "YTGGCNNNNDKCCAR" | "NFIC" | 0.00282 |
| "pos_12" | 2294 | "GCRNRGSMTKCTGGGARWTGTAGTYY" | "Thap11" | 0.000009 |
| "pos_13" | 1276 | "NGCCGCCATCTT" | "YY2" | 0.028003 |
| "pos_14" | 1043 | "NNNNNNNNCAKTTCCTGGAANTNNNNNNNN" | null | null |
| "pos_15" | 749 | "NNVNNNNNRGTCACGTGRCCNNNNNBNNNN" | "TFE3" | 0.001709 |
| "pos_16" | 643 | "TCTCGCGAGA" | "ZBTB33" | 0.000514 |
| "pos_17" | 514 | "GCRGDGSMTKCTGGGAN" | "ZNF76" | 0.000282 |
| "pos_18" | 509 | "NNNNNNNNNNNNNGRTGACMTCATCN" | "JUN" | 0.00138 |
| "pos_19" | 433 | "NCNNNAGGGGGCAGYAGWGAGC" | "CTCF" | 0.002991 |
Most of the common patterns match the motifs of factors known to be active
in K562 cells. The GATA motif of pos_0 is bound there by GATA1, a master
regulator of red blood cell development, and pos_3 and pos_9 are the
MAF recognition element bound by NF-E2, another erythroid factor. Alongside
them are factors found at regulatory elements in most cell types, such as
CTCF, the GC box of the SP and KLF families, ETS factors, NF-Y, NRF1, and
the AP-1 and ATF families. Note that the best match is only the closest
motif in JASPAR, and the factors of a family often share one motif, so
SP4 here stands for the SP and KLF families as a whole rather than SP4
itself. The negative patterns are next:
summary.filter(pl.col("pattern").str.starts_with("neg"))
| pattern | seqlets | consensus | tf | qvalue |
|---|---|---|---|---|
| str | i64 | str | str | f64 |
| "neg_0" | 503 | "CTGATTGGTSCRT" | "NFYA" | 0.039402 |
| "neg_1" | 235 | "CTGTCCRTGGTGCTGAA" | "REST" | 3.3972e-9 |
| "neg_2" | 203 | "NNDNNNDNNWTGATKGATGSATNNATNN" | "PBX1" | 0.03845 |
| "neg_3" | 184 | "VDGNTGGCATCNN" | null | null |
| "neg_4" | 156 | "GYRTGGTGGYGBRCACCTGTART" | "SNAI2" | 0.002955 |
| "neg_5" | 148 | "NAGRTGGGCAYCN" | null | null |
| "neg_6" | 131 | "NNNNATCCATGATGVTNNNNNNN" | null | null |
| "neg_7" | 115 | "NNNNNNATGATGATSATGATGATGATGATG" | null | null |
| "neg_8" | 103 | "GATGAGATTTGGGTGGGGACA" | "KLF4" | 0.018185 |
| "neg_9" | 84 | "GATGAGATTTGGGTGGG" | null | null |
| "neg_10" | 79 | "TNTGATTGGHWNAW" | "NFYA" | 0.002573 |
| "neg_11" | 57 | "KSCRNGNTGGCATSCAYSTG" | "NFIX" | 0.039737 |
| "neg_12" | 54 | "NVNNSATGVNNNNGMTGGCAYCNNNNANAN" | null | null |
| "neg_13" | 54 | "CGGGTTCACGCCATTCTCCTGCCTCAG" | null | null |
| "neg_14" | 49 | "CCACTAGRTGGCAS" | "CTCF" | 0.000026 |
| "neg_15" | 37 | "GCTGAGGCAGGAGAATGGCATGAACC" | null | null |
The negative patterns are much rarer, the most common seen in 503
seqlets, and fewer of them match a known motif. The clearest is neg_1,
the motif of REST, a repressor that silences neuronal genes outside the
nervous system, so the model has learned that REST sites make a region
less accessible in these cells. Others are harder to read: neg_0 and
neg_10 resemble the NF-Y motif, which is also among the positive
patterns, and some, such as the ATGATG repeats of neg_7, match nothing
at all.
Contribution scores along the genome¶
The patterns summarize what the model learned across all regions, while the contribution scores show it at work in each one. The regions file lists the 1 kb windows the scores were computed for:
regions = pl.read_csv(regions_path, separator="\t", has_header=False,
new_columns=["chrom", "start", "end"])
print(regions.height)
regions.head(3)
269718
| chrom | start | end |
|---|---|---|
| str | i64 | i64 |
| "chr1" | 100028707 | 100029707 |
| "chr1" | 100027505 | 100028505 |
| "chr1" | 100028104 | 100029104 |
We look at hypersensitive site 2 of the Ī²-globin locus control region, a well-studied erythroid enhancer that drives the expression of the globin genes in K562 cells, and pick the window that covers its core, keeping its row number in the file for later:
hs2 = regions.with_row_index().filter(
pl.col("chrom") == "chr11", pl.col("start") <= 5_280_600, pl.col("end") >= 5_280_900,
)
hs2
| index | chrom | start | end |
|---|---|---|---|
| u32 | str | i64 | i64 |
| 50341 | "chr11" | 5280325 | 5281325 |
signal_values() reads the score of every
base of the window from the remote track, and hg38() fetches its
sequence. The track holds one score per base, that of the base in the
genome, so view_attributions() draws
the sequence itself with each letter scaled by its contribution, with
positions along the chromosome:
row, chrom, start, end = hs2.row(0)
scores = um.signal_values(scores_url, chrom, start, end)
seq = hg38(chrom, start, end)
print(len(seq), scores.min().round(4), scores.max().round(4))
ax = um.view_attributions(scores, seq, start=start, figsize=(12, 2.5))
1000 -0.0186 0.144
Most of the window contributes little or nothing, apart from a stretch of about 150 bases in the middle that holds the core of HS2, and a few bases score below zero, meaning that the model expects the region to be a little less accessible because of them. To see which patterns the peaks belong to, we scan the window with the positive patterns that matched a known motif and label each hit with its match. A strict P-value cutoff, and removing overlapping hits across strands and patterns, keeps only the clearest sites:
named = {m.name: m for m in ppms}
labels = dict(best.select("subject", "tf").iter_rows())
hits = um.scan_sequences(
[named[p] for p in labels if p.startswith("pos")], seq, pvalue=1e-5,
no_overlaps=True, no_overlaps_by_strand=True, no_overlaps_by_motif=True,
)
hits = hits.with_columns(motif=pl.col("motif").replace_strict(labels)).sort("start")
hits.select("motif", "start", "end", "strand", "match", "pvalue")
| motif | start | end | strand | match | pvalue |
|---|---|---|---|---|---|
| str | i64 | i64 | str | str | f64 |
| "MAF::NFE2" | 426 | 436 | "-" | "GATGAGTCAT" | 0.000003 |
| "MAF::NFE2" | 436 | 446 | "-" | "GCTGAGTCAT" | 9.5367e-7 |
| "SP4" | 501 | 513 | "-" | "TGCCCCACCCCC" | 0.000005 |
The hits are counted along the sequence from 0, as view_attributions()
expects, while region is given in positions along the chromosome, the
same as the axis:
ax = um.view_attributions(scores, seq, start=start, region=(start + 415, start + 520), hits=hits)
The tallest peak is a pair of NF-E2 sites in tandem, GATGAGTCAT and
GCTGAGTCAT on the reverse strand. Together they make up most of the
window's total score:
round(float(scores[426:446].sum() / scores.sum()), 2)
0.66
The tandem NF-E2 sites are the best-known feature of HS2 and are needed for
its activity as an enhancer, and the model has found them from
accessibility alone. The smaller peak to the right is a GC-rich SP or KLF
site, and the one at 5,280,695 in the full window sits on AGATAG, a GATA
site too weak to pass the scan's cutoff.
Hypothetical contribution scores¶
The track only holds the score of the base in the genome at each position.
Attribution methods also give hypothetical scores, the contribution each of
the four letters would make at each position, which show what the model
would accept at a site and not only what is there. TF-MoDISco keeps these
for every seqlet, along with the row of the regions file the seqlet came
from and where it lies in the central 500 bases of that window. Seqlets
found on the reverse strand are stored reverse complemented, so we turn
those back, which for letters in the order ACGT is a flip of both axes.
This reads the seqlets of every pattern that fall in the HS2 window, with
h5py:
import h5py
def seqlets_in(path, row, offset):
found = []
with h5py.File(path, "r") as f:
for group in ("pos_patterns", "neg_patterns"):
for name, pattern in f[group].items():
seqlets = pattern["seqlets"]
for i in np.flatnonzero(seqlets["example_idx"][()] == row):
scores = seqlets["hypothetical_contribs"][i]
if seqlets["is_revcomp"][i]:
scores = scores[::-1, ::-1]
found.append({
"pattern": f"{group[:3]}_{name.rsplit('_', 1)[1]}",
"start": offset + int(seqlets["start"][i]),
"end": offset + int(seqlets["end"][i]),
"scores": scores.T,
})
return found
found = seqlets_in(modisco_path, row, start + 250)
[(s["pattern"], s["start"], s["end"]) for s in found]
[('pos_0', 5280678, 5280708), ('pos_3', 5280751, 5280781)]
Two seqlets were found in the window. One belongs to the NF-E2 pattern and
covers the tandem sites, and the other belongs to the GATA pattern and
covers the AGATAG site that our scan missed, so TF-MoDISco counted it
after all. As a check, the hypothetical scores of the bases in the genome
should be the scores in the bigWig track:
nfe2 = next(s for s in found if s["pattern"] == "pos_3")
a, b = nfe2["start"] - start, nfe2["end"] - start
observed = (um.one_hot(seq[a:b])[0] * nfe2["scores"]).sum(axis=0)
np.allclose(observed, scores[a:b], atol=1e-3)
True
With kind="all", view_attributions() draws every letter at every
position, those that would raise the prediction above zero and those that
would lower it below, and with kind="heatmap" it shows the same scores as
colours, with dots marking the sequence:
fig, axes = plt.subplots(2, 1, figsize=(9, 5.5))
um.view_attributions(nfe2["scores"], seq[a:b], kind="all", start=nfe2["start"],
ylabel="Hypothetical score", ax=axes[0])
um.view_attributions(nfe2["scores"], seq[a:b], kind="heatmap", start=nfe2["start"],
ylabel="Hypothetical score", ax=axes[1])
fig.tight_layout()
Across the two sites, the letter in the genome is nearly always the one the model scores highest, and most of the other letters score below zero, so the model expects almost any change to the sites to make the enhancer less accessible. The scores fade quickly past the second site, from about 5,280,772 on, where no letter matters much.
Motifs as convolution filters¶
The first layer of most sequence models is a convolution, which slides a
set of filters along the sequence much as a scan slides a position weight
matrix. one_hot() gives sequences as arrays of
shape (sequences, letters, length), the layout PyTorch's Conv1d takes,
and to_filters() gives motifs as weights of
shape (filters, letters, width). The pattern's PPM has letters that no
seqlet had, whose log-odds would be infinite, so
normalize() adds a pseudocount first, and
with rc the reverse complement of the motif follows it as a second
filter:
mare = um.normalize(named["pos_3"])
filters = um.to_filters(mare, rc=True)
x = um.one_hot(seq)
print(x.shape, filters.shape)
(1, 4, 1000) (2, 4, 10)
A convolution without a bias sums the product of the filter and the sequence over every window, which numpy can do directly:
windows = np.lib.stride_tricks.sliding_window_view(x, filters.shape[2], axis=-1)
conv = np.einsum("fcw,ncpw->nfp", filters, windows)[0]
print(conv.argmax(axis=1), conv.max(axis=1).round(3))
um.scan_sequences(mare, seq, pvalue=1e-4).select("start", "end", "strand", "score", "pvalue")
[425 436] [10.501 16.835]
| start | end | strand | score | pvalue |
|---|---|---|---|---|
| i64 | i64 | str | f64 | f64 |
| 425 | 435 | "+" | 10.498 | 0.000074 |
| 426 | 436 | "-" | 15.941 | 0.000003 |
| 436 | 446 | "-" | 16.83 | 9.5367e-7 |
The best window of the forward filter starts at 425 and that of the
reverse filter at 436, the same sites that
scan_sequences() finds, with the same
scores to within a few thousandths, since scanning rounds the matrix to
compute exact P-values. The same filters can be given to a model, for
example to start its first layer from known motifs:
import torch
conv = torch.nn.Conv1d(4, filters.shape[0], filters.shape[2], bias=False)
with torch.no_grad():
conv.weight.copy_(torch.from_numpy(filters))
Going the other way, from_filters() turns
filters into motifs, so the first-layer weights of a trained model can be
compared with known motifs or drawn like any other. Here we turn our two
filters back into motifs:
back = um.from_filters(filters, names=["MARE", "MARE, reverse complement"])
fig = um.view_motifs(back, rc=False, min_overlap=100, names_pos="right", figsize=(6, 3.5))
Reading contribution score files¶
Here the scores came from a bigWig track and from TF-MoDISco's seqlets, but
the full hypothetical scores of every region are usually kept in files of
their own. chromBPNet writes them, along with the one-hot sequences, to
HDF5 files named <prefix>.counts_scores.h5 and
<prefix>.profile_scores.h5, next to a file of the regions they cover,
and for a full set of peaks these are often tens of gigabytes.
read_contributions() reads them,
optionally just some of the sequences, and gives each sequence its genomic
window, so that its view() method labels positions along the chromosome:
contribs = um.read_contributions("model.counts_scores.h5", indices=range(100))
contribs.names[:3]
contribs.view(0, kind="heatmap")
The scores and one-hot sequences that TF-MoDISco itself takes, as numpy files, are read the same way:
contribs = um.read_contributions("shap.npz", sequences="ohe.npz")
From there, the scores of each sequence are a numpy array, and the sequences are strings that can be scanned, compared or summarized with the rest of the package.