Scanning and enrichment¶
This guide covers finding motif sites in sequences and asking which motifs stand out in them: scanning, choosing a P-value cutoff, dealing with both strands and with overlapping hits, testing motifs for enrichment, for co-occurrence and for positional preferences, and finally discovering motifs de novo. Preparing the sequences and backgrounds that go into these analyses is covered in Working with sequences.
The example data was made with these analyses in mind. There are twelve example motifs, a hundred AT-rich target sequences of 1000 bases each, and a hundred control sequences just like them. W-box sites were planted in 60 of the targets and G-box sites in 40, towards the middle of the sequences, while nothing was planted in the controls. Throughout this guide, the question is how well the W-box and the G-box can be picked out.
import warnings
import polars as pl
import universalmotif as um
motifs = um.example_motifs()
targets = um.example_sequences("targets")
controls = um.example_sequences("controls")
print([m.name for m in motifs])
['TATA-box', 'AP-1', 'CRE', 'G-box', 'E-box', 'W-box', 'GATA', 'GC-box', 'ETS', 'CArG-box', 'Telo-box', 'CCAAT-box']
Scanning sequences¶
scan_sequences() scores every window the
width of the motif in every sequence, and reports those whose P-value is at
or below a cutoff, which is 1e-4 by default. Let us scan the targets for the
W-box:
by_name = {m.name: m for m in motifs}
wbox = by_name["W-box"]
hits = um.scan_sequences(wbox, targets)
print(hits.height)
hits.select("sequence", "start", "end", "strand", "score", "score_pct", "match", "pvalue").head()
30
| sequence | start | end | strand | score | score_pct | match | pvalue |
|---|---|---|---|---|---|---|---|
| str | i64 | i64 | str | f64 | f64 | str | f64 |
| "target_004" | 448 | 457 | "-" | 12.341 | 89.832535 | "GTTTGACCG" | 0.000034 |
| "target_005" | 275 | 284 | "+" | 12.341 | 89.832535 | "GTTTGACCG" | 0.000034 |
| "target_008" | 646 | 655 | "-" | 12.476 | 90.815226 | "GGTTGACCG" | 0.000031 |
| "target_015" | 348 | 357 | "-" | 10.95 | 79.707176 | "CGTTGACTC" | 0.000099 |
| "target_017" | 292 | 301 | "+" | 12.341 | 89.832535 | "GTTTGACCG" | 0.000034 |
Each row is one hit. start and end count from 0 and leave the end out,
as Python slices do, so a hit's letters are sequence[start:end]. They are
always given on the forward strand, while match shows the site in the
orientation it matched in, so hits on the - strand are reverse
complemented. score_pct gives the score as a percentage of the highest
score the motif can reach. Not shown here are the motif_i and
sequence_i columns, which are 0-based indices into the motifs and
sequences that were scanned.
Choosing a P-value cutoff¶
The score of a site is the sum of the motif's log-odds weights for its
letters, and its P-value is the probability that a random site would score
at least as much. Each motif's P-value cutoff is thus turned into a score
threshold, which motif_pvalue() can
calculate directly. For the W-box at three different cutoffs:
print(um.motif_pvalue(wbox, pvalue=[1e-3, 1e-4, 1e-5]))
print(um.motif_range(wbox))
[ 6.845 10.885 13.306] (-37.998, 13.733)
A site must score at least 10.885 to pass the default cutoff of 1e-4, out of a possible range of -37.998 to 13.733. At 1e-5 the threshold rises to 13.306, close to the maximum, so only near-perfect sites would pass.
So which cutoff to use? It helps to consider how many hits to expect by chance. Scanning a 9-base motif across 100 sequences of 1000 bases on both strands means scoring 100 Ć 992 Ć 2 = 198,400 windows, so a cutoff of 1e-4 should let through about 20 hits in purely random sequence, and a cutoff of 1e-3 about 200. Let us compare the hits in the targets and the controls at a few cutoffs:
rows = []
for cutoff in (1e-3, 1e-4, 1e-5):
t = um.scan_sequences(wbox, targets, pvalue=cutoff)
c = um.scan_sequences(wbox, controls, pvalue=cutoff)
rows.append({
"cutoff": cutoff,
"target_hits": t.height, "target_seqs": t["sequence"].n_unique(),
"control_hits": c.height, "control_seqs": c["sequence"].n_unique(),
})
pl.DataFrame(rows)
| cutoff | target_hits | target_seqs | control_hits | control_seqs |
|---|---|---|---|---|
| f64 | i64 | i64 | i64 | i64 |
| 0.001 | 153 | 81 | 93 | 65 |
| 0.0001 | 30 | 30 | 10 | 9 |
| 0.00001 | 1 | 1 | 0 | 0 |
At 1e-3, most sequences in both sets have a hit, 81 of the targets and 65 of the controls, so the planted sites are drowned out by chance matches. At 1e-4 the difference is much clearer, with 30 target sequences against 9 controls, though this also means that only half of the 60 targets with a planted W-box have a hit at all. The planted sites were drawn from the motif's probabilities, so many of them are imperfect, and a strict cutoff misses them. At 1e-5 almost nothing is left.
It is important to remember that there is likely never one single cutoff which works in every situation. It depends on the motif, on the sequences, and on whether missing real sites or reporting false ones is the bigger problem for what you are doing. 1e-4 is a reasonable place to start for motifs of this size, while longer and more informative motifs can afford stricter cutoffs. Note also that the controls have only 10 hits at 1e-4, half of the 20 we expected; we will come back to this shortly.
Short motifs have a further limit: there are only so many possible sites, so the best possible P-value is not that small. Let us calculate it for each motif from its highest possible score:
best = [um.motif_pvalue(m, score=um.motif_range(m)[1]) for m in motifs]
with pl.Config(tbl_rows=12):
display(pl.DataFrame({
"motif": [m.name for m in motifs],
"width": [m.width for m in motifs],
"best_pvalue": best,
}))
| motif | width | best_pvalue |
|---|---|---|
| str | i64 | f64 |
| "TATA-box" | 10 | 9.5367e-7 |
| "AP-1" | 9 | 0.000008 |
| "CRE" | 8 | 0.000015 |
| "G-box" | 8 | 0.000015 |
| "E-box" | 8 | 0.000015 |
| "W-box" | 9 | 0.000004 |
| "GATA" | 8 | 0.000031 |
| "GC-box" | 9 | 0.000004 |
| "ETS" | 10 | 9.5367e-7 |
| "CArG-box" | 10 | 9.5367e-7 |
| "Telo-box" | 9 | 0.000004 |
| "CCAAT-box" | 9 | 0.000004 |
The best P-value of a DNA motif of width w can be no smaller than 4 to the power of -w, and here they range from about 1e-6 to 3e-5. The shortest motifs, such as the G-box, cannot reach a cutoff of 1e-5 at all, and scanning with one returns no hits for them, with a warning saying so:
with warnings.catch_warnings(record=True) as caught:
print(um.scan_sequences(by_name["G-box"], targets, pvalue=1e-5).height)
print(caught[0].message)
0 1 motif cannot reach a P-value of 1e-05 and will have no hits; the smallest P-value each can give is: G-box (1.53e-05)
The background matters¶
The P-values above are calculated against each motif's background, which for the example motifs is uniform:
wbox.bkg
{'A': 0.25, 'C': 0.25, 'G': 0.25, 'T': 0.25}
The example sequences however are AT-rich, so motifs made mostly of As
and Ts will match them more often than their P-values suggest, and motifs
rich in Cs and Gs less often. The fix is to give the motifs the background
of the sequences being scanned, which we can calculate with
get_bkg(). Let us compare the number of hits in
the targets for every motif before and after:
bkg = um.get_bkg(targets, k=1)
bkg = dict(zip(bkg["kmer"], bkg["probability"]))
matched = [m.replace(bkg=bkg) for m in motifs]
def hits_per_motif(ms, seqs):
return um.scan_sequences(ms, seqs).group_by("motif").len()
with pl.Config(tbl_rows=12):
display(
hits_per_motif(motifs, targets).rename({"len": "uniform"})
.join(hits_per_motif(matched, targets).rename({"len": "matched"}), on="motif")
.sort("motif")
)
| motif | uniform | matched |
|---|---|---|
| str | u32 | u32 |
| "AP-1" | 24 | 25 |
| "CArG-box" | 24 | 18 |
| "CCAAT-box" | 12 | 17 |
| "CRE" | 6 | 14 |
| "E-box" | 5 | 16 |
| "ETS" | 8 | 20 |
| "G-box" | 27 | 59 |
| "GATA" | 30 | 26 |
| "GC-box" | 3 | 26 |
| "TATA-box" | 111 | 11 |
| "Telo-box" | 36 | 22 |
| "W-box" | 30 | 47 |
The TATA-box, made mostly of As and Ts, goes from 111 hits down to 11, while GC-richer motifs such as the GC-box and the G-box gain hits. The 111 TATA-box hits were more than five times the 20 or so expected by chance, and the W-box had half as many hits in the controls as expected; both were symptoms of the mismatched background. We will use the motifs with matched backgrounds from here on.
Both strands¶
DNA motifs can bind either strand, so by default both are scanned. A palindromic motif such as the G-box, CACGTG, reads the same on both strands, which means each of its sites is found twice:
by_name = {m.name: m for m in matched}
wbox, gbox = by_name["W-box"], by_name["G-box"]
gbox_hits = um.scan_sequences(gbox, targets)
print(gbox_hits.height)
gbox_hits.select("sequence", "start", "end", "strand", "match", "pvalue").head(6)
59
| sequence | start | end | strand | match | pvalue |
|---|---|---|---|---|---|
| str | i64 | i64 | str | str | f64 |
| "target_004" | 716 | 724 | "+" | "GCACGTGT" | 0.000043 |
| "target_004" | 716 | 724 | "-" | "ACACGTGC" | 0.000079 |
| "target_008" | 471 | 479 | "+" | "GCACGTGT" | 0.000043 |
| "target_008" | 471 | 479 | "-" | "ACACGTGC" | 0.000079 |
| "target_015" | 187 | 195 | "+" | "TCACGTGA" | 0.000038 |
| "target_015" | 187 | 195 | "-" | "TCACGTGA" | 0.000038 |
Each G-box site turns up once per strand at the same position. The two
P-values can differ a little, since the flanks of the motif are not quite
symmetrical, but they are the same site. Scanning only the forward strand
with rc=False gives about half as many hits:
um.scan_sequences(gbox, targets, rc=False).height
31
Of course, for motifs which are not palindromic, scanning only the forward strand would miss the sites on the other.
Overlapping hits¶
Hits can also overlap one another, either because of palindromes as above
or because a sequence has several good matches in a row, as happens with
low-complexity sequence. dedup_hits() keeps
only the best of each group of overlapping hits, by P-value by default. Hits
on different strands are kept apart unless ignore_strand=True, and hits of
different motifs unless ignore_motif=True:
print(um.dedup_hits(gbox_hits).height)
print(um.dedup_hits(gbox_hits, ignore_strand=True).height)
58 33
Deduplicating within strands removes a single hit, whereas ignoring the strands brings the 59 hits down to 33, one per site. Overlapping hits on the same strand are more of a problem for motifs which happen to suit the composition of the sequences. Let us compare the TATA-box with its uniform and its matched background, printing the number of hits before deduplication, after deduplicating within strands, and after deduplicating across them:
tata = by_name["TATA-box"]
tata_uniform = motifs[0]
for m in (tata_uniform, tata):
h = um.scan_sequences(m, targets)
print(h.height, um.dedup_hits(h).height, um.dedup_hits(h, ignore_strand=True).height)
111 104 94 11 11 11
With the uniform background, the TATA-box finds 111 hits in the AT-rich
targets, and 17 of them overlap better hits on the same strand or the
other. With the matched background, none of its 11 hits overlap. Setting
no_overlaps=True in scan_sequences()
deduplicates as part of the scan.
Enrichment against control sequences¶
Now, to the question of which motifs are enriched in the targets.
enrich_motifs() scans the targets and the
background sequences, then tests each motif with a one-sided Fisher's exact
test. By default the test is on the number of sequences with at least one
hit, and only motifs with a Benjamini-Hochberg q-value of at most 0.1 are
returned:
cols = ["motif", "target_seq_hits", "bkg_seq_hits", "target_site_hits",
"bkg_site_hits", "enrichment", "pvalue", "qvalue"]
um.enrich_motifs(matched, targets, controls).select(cols)
| motif | target_seq_hits | bkg_seq_hits | target_site_hits | bkg_site_hits | enrichment | pvalue | qvalue |
|---|---|---|---|---|---|---|---|
| str | i64 | i64 | i64 | i64 | f64 | f64 | f64 |
| "G-box" | 32 | 7 | 59 | 10 | 4.571429 | 0.000009 | 0.000103 |
| "W-box" | 46 | 19 | 47 | 21 | 2.421053 | 0.000046 | 0.000275 |
Both planted motifs come out on top, and nothing else passes. The G-box has
a hit in 32 targets against 7 controls, and the W-box in 46 against 19. The
numbers of hits are given too, but note that every site of a palindromic
motif such as the G-box counts twice. Setting qvalue=1 returns every
motif, which shows how the rest fared:
um.enrich_motifs(matched, targets, controls, qvalue=1).select(cols).head(5)
| motif | target_seq_hits | bkg_seq_hits | target_site_hits | bkg_site_hits | enrichment | pvalue | qvalue |
|---|---|---|---|---|---|---|---|
| str | i64 | i64 | i64 | i64 | f64 | f64 | f64 |
| "G-box" | 32 | 7 | 59 | 10 | 4.571429 | 0.000009 | 0.000103 |
| "W-box" | 46 | 19 | 47 | 21 | 2.421053 | 0.000046 | 0.000275 |
| "GC-box" | 21 | 12 | 26 | 13 | 1.75 | 0.068299 | 0.273196 |
| "AP-1" | 14 | 10 | 25 | 17 | 1.4 | 0.264734 | 0.698521 |
| "Telo-box" | 19 | 15 | 22 | 16 | 1.266667 | 0.29105 | 0.698521 |
Note that the W-box now has 21 hits in the controls, close to the 20 we
expected by chance, since the background matches. Here too the P-value
cutoff matters. At a cutoff of 1e-3, so many sequences in both sets have a
hit that the per-sequence test finds nothing at all, while counting sites
with test="sites" still picks out the planted motifs, along with the CRE
(which shares a TGAC with the W-box and an ACGT with the G-box):
print(um.enrich_motifs(matched, targets, controls, pvalue=1e-3).height)
um.enrich_motifs(matched, targets, controls, pvalue=1e-3, test="sites").select(cols)
0
| motif | target_seq_hits | bkg_seq_hits | target_site_hits | bkg_site_hits | enrichment | pvalue | qvalue |
|---|---|---|---|---|---|---|---|
| str | i64 | i64 | i64 | i64 | f64 | f64 | f64 |
| "W-box" | 94 | 88 | 258 | 185 | 1.394595 | 0.000311 | 0.003731 |
| "CRE" | 82 | 74 | 236 | 172 | 1.372093 | 0.000909 | 0.005453 |
| "G-box" | 86 | 79 | 246 | 193 | 1.274611 | 0.00657 | 0.026278 |
The site test compares the number of hits out of every position a motif could occupy, so it stays sensitive when most sequences have a hit, but it is also more easily swayed by a few sequences with many hits, such as repeats. Which test suits your data better is a judgment call; for short sequences such as ChIP-seq peaks, where a site is generally either there or not, the per-sequence test is usually the more natural choice.
Enrichment against shuffled sequences¶
Matching control sequences are not always available. Without
bkg_sequences, enrich_motifs() makes its
own background by shuffling the targets with
shuffle_sequences(), keeping their
dinucleotide composition by default:
um.enrich_motifs(matched, targets, seed=1).select(cols)
| motif | target_seq_hits | bkg_seq_hits | target_site_hits | bkg_site_hits | enrichment | pvalue | qvalue |
|---|---|---|---|---|---|---|---|
| str | i64 | i64 | i64 | i64 | f64 | f64 | f64 |
| "W-box" | 46 | 18 | 47 | 19 | 2.555556 | 0.000022 | 0.000267 |
| "G-box" | 32 | 16 | 59 | 30 | 2.0 | 0.007069 | 0.042415 |
Both planted motifs are found again, though the G-box less convincingly
than against the controls, since 16 of the shuffled sequences have a hit
against 7 of the controls. A shuffled background keeps the composition of
the targets, which makes it a reasonable null when nothing better is at
hand, but real control sequences share far more with the targets and are
the better choice whenever they exist. The results also depend somewhat on
the seed and on shuffle_k, so it is worth checking that a result does
not hinge on one particular shuffle.
Co-occurrence¶
Transcription factors often work together, so beyond asking whether each
motif is enriched on its own, one may ask whether pairs of motifs turn up
in the same sequences more often than chance would predict.
motif_coocc() scans the sequences, counts
for every pair of motifs the sequences with hits from one, the other, both
or neither, and runs a one-sided Fisher's exact test on each:
coocc = um.motif_coocc(matched, targets)
coocc.select("motif_a", "motif_b", "a_only", "b_only", "both", "neither", "pvalue", "qvalue").head(5)
| motif_a | motif_b | a_only | b_only | both | neither | pvalue | qvalue |
|---|---|---|---|---|---|---|---|
| str | str | i64 | i64 | i64 | i64 | f64 | f64 |
| "W-box" | "CArG-box" | 35 | 3 | 11 | 51 | 0.008879 | 0.550477 |
| "AP-1" | "W-box" | 4 | 36 | 10 | 50 | 0.038054 | 0.898635 |
| "GC-box" | "ETS" | 14 | 11 | 7 | 68 | 0.046138 | 0.898635 |
| "ETS" | "CArG-box" | 13 | 9 | 5 | 73 | 0.074914 | 0.898635 |
| "TATA-box" | "GATA" | 6 | 15 | 4 | 75 | 0.09275 | 0.898635 |
None of the pairs passes correction, which is reassuring since the example
data was not made with co-occurrence in mind. To see what a real
co-occurrence looks like, let us make 100 new AT-rich sequences, plant an
AP-1 site at position 200 and an E-box site 25 bases further along in the
first 40 of them, and test AP-1, the E-box and GATA. Setting
max_distance adds two descriptive columns: the number of shared
sequences in which the two motifs have hits starting within that distance
of each other, and the median of the nearest distances:
hosts = um.create_sequences(seqnum=100, seqlen=500, freqs=bkg, seed=3)
in_first_40 = [[200] if i < 40 else [] for i in range(100)]
paired = um.implant_motifs(by_name["AP-1"], hosts, positions=in_first_40, seed=4)
paired = um.implant_motifs(by_name["E-box"], paired,
positions=[[p + 25 for p in ps] for ps in in_first_40], seed=5)
trio = [by_name["AP-1"], by_name["E-box"], by_name["GATA"]]
um.motif_coocc(trio, paired, max_distance=50).select(
"motif_a", "motif_b", "both", "pvalue", "qvalue", "both_clustered", "median_distance"
)
| motif_a | motif_b | both | pvalue | qvalue | both_clustered | median_distance |
|---|---|---|---|---|---|---|
| str | str | i64 | f64 | f64 | i64 | f64 |
| "AP-1" | "E-box" | 11 | 0.004205 | 0.012615 | 10 | 25.0 |
| "AP-1" | "GATA" | 2 | 0.779078 | 0.87809 | 2 | 21.5 |
| "E-box" | "GATA" | 1 | 0.87809 | 0.87809 | 1 | 4.0 |
AP-1 and the E-box come out as co-occurring, and in 10 of the 11 sequences
with both, their hits are within 50 bases of each other at a median
distance of 25, just as planted. The pair is only found in 11 of the 40
sequences it was planted in though, which is once more a matter of the
cutoff, as many planted sites are too degenerate to reach a P-value of
1e-4. The E-box and GATA share no sequences at all, so with the default
min_coocc=1 they are not tested and their P-value is missing. The
distance columns do not change the test itself, and are best thought of as
an aid to interpretation, for example to spot pairs of motifs which might
bind together as a complex.
Positional enrichment¶
Sites of a motif may also favour a particular region of the sequences, such
as the middle of ChIP-seq peaks or the region just upstream of
transcription start sites. motif_peaks()
takes a table of hits from sequences of equal length and tests, for each
motif, whether the centres of its best hit per sequence fall in a window
around the middle of the sequences more often than they would if they were
spread evenly. Several window sizes are tried, and the best is reported:
hits = um.scan_sequences(matched, targets)
peaks = um.motif_peaks(hits, seq_length=1000)
peaks.select("motif", "nhits", "best_window", "best_centre", "hits_in",
"expected_in", "pvalue", "qvalue")
| motif | nhits | best_window | best_centre | hits_in | expected_in | pvalue | qvalue |
|---|---|---|---|---|---|---|---|
| str | i64 | i64 | f64 | i64 | f64 | f64 | f64 |
| "W-box" | 46 | 490 | 499.0 | 37 | 22.54 | 0.000568 | 0.006811 |
Only the W-box passes. Its best window is 490 bases wide and centred on the
middle of the sequences, holding 37 of the 46 best hits against the 22.5
expected for a window of that size. Such a wide window suggests the
planted sites are spread across a broad middle region rather than tightly
centred. Let us plot the hit positions of the W-box and the G-box with
plot_motif_peaks(), setting
qvalue=1 to keep every motif in the results:
peaks = um.motif_peaks(hits, seq_length=1000, qvalue=1)
display(peaks.filter(pl.col("motif").is_in(["W-box", "G-box"]))
.select("motif", "nhits", "best_window", "hits_in", "pvalue", "qvalue"))
fig = um.plot_motif_peaks(peaks, motifs=["W-box", "G-box"], figsize=(8, 3))
| motif | nhits | best_window | hits_in | pvalue | qvalue |
|---|---|---|---|---|---|
| str | i64 | i64 | i64 | f64 | f64 |
| "W-box" | 46 | 490 | 37 | 0.000568 | 0.006811 |
| "G-box" | 32 | 480 | 25 | 0.024327 | 0.145961 |
The G-box hits lean towards the middle as well, with 25 of its 32 best hits
in a window of 480 bases, but with fewer sequences carrying a hit the
evidence is weaker, and its q-value of 0.146 does not pass the default
cutoff of 0.1. For sequences without a natural centre, mode="local" moves
the window along the whole length of the sequences instead of keeping it in
the middle.
Motif match matrices¶
For many sequences at once, such as the peaks of a single-cell ATAC
experiment, it is often enough to know which motifs match which sequences.
match_motifs() scans as
scan_sequences() does, but keeps only a
sparse matrix with a row per sequence and a column per motif, holding a 1
for each match, the number of hits with value="counts", or the best score
with value="scores". Since the hits themselves are never kept, this works
for hundreds of thousands of peaks. It needs the scipy extra:
matches = um.match_motifs(motifs, targets, value="counts")
print(matches)
pl.DataFrame({
"motif": matches.motifs,
"sequences": (matches.matrix > 0).sum(axis=0),
"hits": matches.matrix.sum(axis=0),
}).sort("sequences", descending=True).head(4)
MotifMatches(100 sequences x 12 motifs, 120 counts)
| motif | sequences | hits |
|---|---|---|
| str | i64 | u64 |
| "TATA-box" | 38 | 55 |
| "Telo-box" | 18 | 23 |
| "W-box" | 17 | 17 |
| "G-box" | 12 | 18 |
The P-value cutoff defaults to 5e-5 here, stricter than the 1e-4 of
scan_sequences(), as is usual for
motif annotations of peaks. With
to_anndata() the matrix becomes an AnnData object, and
add_motif_matches() stores it in an
existing AnnData object of peaks in the form that pychromVAR's
compute_deviations() reads, for chromVAR's motif accessibility scores.
De novo discovery¶
So far we have known which motifs to look for. Often though the question is
which motifs are present at all, and that calls for de novo discovery.
motif_finder() looks, at each motif width
in a range, for short words found in more target than background sequences,
builds motifs around the best of them and refines them, keeping those whose
hits are enriched and then masking their sites to look for more. Let us run
it on the targets against the controls, keeping the widths to a modest
range:
found = um.motif_finder(targets, controls, min_width=6, max_width=10, seed=1, nthreads=2)
found.select("name", "consensus", "width", "seqs_pos", "seqs_neg", "candidates", "pvalue", "qvalue")
| name | consensus | width | seqs_pos | seqs_neg | candidates | pvalue | qvalue |
|---|---|---|---|---|---|---|---|
| str | str | i64 | i64 | i64 | i64 | f64 | f64 |
| "motif_1" | "CACGTG" | 6 | 33 | 4 | 7994 | 0.000351 | 0.014892 |
Only one motif passes the default q-value cutoff of 0.05, and it is clearly
the G-box, found in 33 targets and only 4 controls. Its P-value comes from
Fisher's exact test on the sequences with hits, multiplied by the number of
candidate words tried in the search that found it, given in candidates,
since a motif picked as the best of thousands of words would otherwise look
far more significant than it is. The q-values likewise count every word
tried at every width. Setting qvalue=1 shows the candidates that did not
make it:
everything = um.motif_finder(targets, controls, min_width=6, max_width=10, qvalue=1, seed=1, nthreads=2)
everything.select("name", "consensus", "seqs_pos", "seqs_neg", "pvalue", "qvalue")
| name | consensus | seqs_pos | seqs_neg | pvalue | qvalue |
|---|---|---|---|---|---|
| str | str | i64 | i64 | f64 | f64 |
| "motif_1" | "CACGTG" | 33 | 4 | 0.000351 | 0.014892 |
| "motif_2" | "CGACAG" | 33 | 11 | 0.284159 | 1.0 |
| "motif_3" | "AGTGCG" | 24 | 6 | 0.590302 | 1.0 |
| "motif_4" | "GACGTAC" | 16 | 1 | 0.615604 | 1.0 |
| "motif_5" | "AMGACAG" | 31 | 10 | 1.0 | 1.0 |
| "motif_6" | "AGGAARAA" | 29 | 10 | 1.0 | 1.0 |
| "motif_7" | "ATTAAWTTC" | 30 | 11 | 1.0 | 1.0 |
| "motif_8" | "AAAATACWT" | 30 | 11 | 1.0 | 1.0 |
None of the others comes close. The second, in 33 targets against 11
controls, would look convincing on its own, but among thousands of words
some will show a contrast like that by chance. Comparing all of them with
the example motifs using
compare_motifs(), keeping the best match
for each, confirms that only the first resembles anything that was planted:
new_motifs = list(everything["motif"])
n = len(new_motifs)
comparison = um.compare_motifs(new_motifs + motifs, compare_to=list(range(n)), qvalue=1)
(
comparison.filter(pl.col("target_i") >= n)
.sort("pvalue")
.group_by("subject", maintain_order=True).first()
.sort("subject")
.select("subject", "subject_consensus", "target", "target_consensus", "score", "pvalue")
)
| subject | subject_consensus | target | target_consensus | score | pvalue |
|---|---|---|---|---|---|
| str | str | str | str | f64 | f64 |
| "motif_1" | "CACGTG" | "G-box" | "CACGTG" | 0.989654 | 0.000787 |
| "motif_2" | "CGACA" | "CRE" | "TGACG" | 0.732412 | 0.035968 |
| "motif_3" | "AGTGC" | "G-box" | "NCACG" | 0.649335 | 0.124642 |
| "motif_4" | "GACGTAC" | "CRE" | "GACGTCA" | 0.645447 | 0.0332 |
| "motif_5" | "AMGACA" | "CRE" | "ACGTCA" | 0.642737 | 0.127195 |
| "motif_6" | "GGAARAA" | "CArG-box" | "WTWWTGG" | 0.446466 | 0.26662 |
| "motif_7" | "ATTAAWTT" | "TATA-box" | "ATAWAWRN" | 0.526546 | 0.193748 |
| "motif_8" | "AAAATACW" | "CArG-box" | "CCWWTWWT" | 0.501572 | 0.204752 |
Even corrected, the P-values remain somewhat optimistic, since each motif was also refined on the very sequences it is tested on. A good way to check is to test the motifs on sequences they have not seen. Let us make a fresh set by planting G-box and W-box sites in 100 new AT-rich sequences, keeping another 100 as controls:
seqs = um.create_sequences(seqnum=200, seqlen=1000, freqs=bkg, seed=10)
new_targets = um.implant_motifs([gbox, wbox], seqs[:100], centre_bias=3, seed=11)
new_controls = seqs[100:]
um.enrich_motifs(new_motifs, new_targets, new_controls, qvalue=1).select(
"motif", "consensus", "target_seq_hits", "bkg_seq_hits", "pvalue", "qvalue"
)
| motif | consensus | target_seq_hits | bkg_seq_hits | pvalue | qvalue |
|---|---|---|---|---|---|
| str | str | i64 | i64 | f64 | f64 |
| "motif_1" | "CACGTG" | 25 | 2 | 0.000002 | 0.000014 |
| "motif_2" | "CGACAG" | 17 | 10 | 0.114302 | 0.457207 |
| "motif_7" | "ATTAAWTTC" | 23 | 20 | 0.367928 | 0.666667 |
| "motif_5" | "AMGACAG" | 13 | 11 | 0.417058 | 0.666667 |
| "motif_3" | "AGTGCG" | 16 | 14 | 0.423811 | 0.666667 |
| "motif_6" | "AGGAARAA" | 20 | 19 | 0.5 | 0.666667 |
| "motif_4" | "GACGTAC" | 8 | 8 | 0.597107 | 0.682408 |
| "motif_8" | "AAAATACWT" | 13 | 15 | 0.723145 | 0.723145 |
Only the G-box holds up, and none of the other candidates is significant in the new sequences; they were chance patterns of the original ones. Comparing the discovered G-box with the original shows how closely it was recovered:
fig = um.view_motifs([new_motifs[0], by_name["G-box"]])
The discovered motif is sharper than the original, likely because it was built from the planted sites which matched best, but its core is the same.
What about the W-box? It was planted in more sequences than the G-box, yet it was not recovered. The reason becomes clearer when counting the sequences which contain a core word of each motif on either strand:
def with_word(seqs, word):
rc = word[::-1].translate(str.maketrans("ACGT", "TGCA"))
return sum(word in s or rc in s for s in seqs.values())
for word in ("CACGTG", "TTGACC"):
print(word, with_word(targets, word), with_word(controls, word))
CACGTG 30 3 TTGACC 48 24
The G-box core CACGTG is in 30 targets and only 3 controls, a contrast which is hard to miss. A core W-box word, TTGACC, is in 48 targets but also in 24 controls, since the W-box sites were planted with more variation and AT-rich sequences hold plenty of similar words by chance. That contrast is no better than those of chance words such as the second candidate found above, in 33 targets against 11 controls, so it cannot stand out during discovery, even though the W-box is clearly enriched once we know to look for it. Degenerate motifs like this one simply need more sequences, or a more focused set of them, before they can be found de novo.
Where next¶
The API reference describes every option of the functions used here, and the Working with sequences guide covers generating, shuffling and masking the sequences that go into them.