Comparing motifs and P-values¶
Two topics come up again and again when working with motifs, though they are
rarely discussed in much detail: how to tell whether a match to a motif is any
good, and how to tell whether two motifs are alike. This guide covers both. The
first half explains motif scores and their P-values, which are what
scan_sequences() uses to decide what counts
as a hit, and the second half covers comparing motifs with one another and
using those comparisons to merge similar motifs and to cluster them into trees.
For a quicker tour of the package as a whole, see the
getting started guide, and for every argument of every
function the API reference.
import warnings
import universalmotif as um
Most of the examples use the twelve example motifs that come with the package. They are synthetic, each built around a well-known transcription factor binding site, and come as count matrices with the family of the factor set. Let us keep them in a dictionary as well, so that they can be picked out by name. We will also need numpy and polars, and let polars show tables of up to 20 rows in full:
import numpy as np
import polars as pl
pl.Config.set_tbl_rows(20)
motifs = um.example_motifs()
by_name = {m.name: m for m in motifs}
um.summarize_motifs(motifs).select("name", "family", "consensus", "icscore")
| name | family | consensus | icscore |
|---|---|---|---|
| str | str | str | f64 |
| "TATA-box" | "TBP" | "KTATAWAWRN" | 8.82611 |
| "AP-1" | "bZIP" | "NTGASTCAN" | 9.646674 |
| "CRE" | "bZIP" | "TGACGTCA" | 10.72504 |
| "G-box" | "bZIP" | "NCACGTGN" | 7.899263 |
| "E-box" | "bHLH" | "NCAGCTGN" | 8.214013 |
| "W-box" | "WRKY" | "GBTTGACYS" | 9.797764 |
| "GATA" | "GATA" | "CWGATART" | 6.740901 |
| "GC-box" | "C2H2 zinc finger" | "GGGGCGGGG" | 11.75613 |
| "ETS" | "ETS" | "NCCGGAAGTG" | 10.570188 |
| "CArG-box" | "MADS" | "CCWWTWWTGG" | 8.986613 |
| "Telo-box" | "Myb-related" | "AAACCCTAA" | 11.240669 |
| "CCAAT-box" | "NF-Y" | "CNCCAATNH" | 7.712225 |
Scoring a match¶
A motif scores a sequence with its position weight matrix, or PWM, in which each entry is the log2 ratio of the probability of a letter at that position to its probability in the background. Let us look at the PWM of the AP-1 motif:
ap1 = by_name["AP-1"]
um.convert_type(ap1, "PWM").matrix.round(2)
array([[-0.29, -1.94, -5.64, 1.69, -5.64, -2.47, -3.32, 1.95, 0.38],
[-0.15, -3.32, -2.47, -1.25, 0.96, -2.47, 1.95, -3.32, -1.56],
[-0.43, -5.64, 1.86, -2.47, 0.96, -2.47, -5.64, -5.64, 0.7 ],
[ 0.62, 1.86, -2.47, -2.47, -3.32, 1.79, -5.64, -5.64, -0.43]])
The score of a word as long as the motif is the sum of the PWM entries for its
letters, one from each column, and score_match()
does that sum for us. For example, here are the scores of the consensus site,
of the same site with its central C changed to a G, of a CRE site, and of a run
of As:
words = ["ATGACTCAT", "ATGAGTCAT", "CTGACGTCA", "AAAAAAAAA"]
um.score_match(ap1, words)
array([ 11.319, 11.319, -4.858, -15.296])
The first two score the same, since the motif has no preference between C and
G at its centre, while the CRE site, which shares the TGAC half of the site but
not the rest, scores well below zero. The range of possible scores comes from
motif_range(), which adds up the lowest and the
highest entry of every column:
um.motif_range(ap1)
(-35.151, 13.356)
Every score is computed with the PWM entries truncated to three decimals, the
same way scanning computes them, which is why the scores so far have at most
three decimal places. A raw score is hard to interpret on its own however. A score of 11.3 is
close to the best this motif can do, but a longer or more informative motif
would have a different range altogether. One way around this is to give
thresholds as fractions of the range of scores, which is what
motif_score() does; here are the scores at 0%,
50%, 80% and 100% of the range:
um.motif_score(ap1, [0, 0.5, 0.8, 1])
array([-35.151 , -10.8975, 3.6546, 13.356 ])
As we will see shortly though, a fraction of the range means rather different things for different motifs. P-values are a much better yardstick.
From scores to P-values¶
The P-value of a score is the probability that a random word, drawn one letter
at a time from the background, scores at least that much. The background stored
in the example motifs is uniform, so every 9-letter word has the same
probability, as prob_match() shows:
um.prob_match(ap1, ["ATGACTCAT", "AAAAAAAAA"])
array([3.81469727e-06, 3.81469727e-06])
This is simply 1/49. For a motif this short we could therefore work
out a P-value by brute force. get_scores() lists
the scores of all 49 = 262,144 possible words, and since they are
all equally likely, the P-value of a score is the fraction of them that reach
it. Let us do this for the score of the consensus site:
all_scores = um.get_scores(ap1)
print(len(all_scores), (all_scores >= 11.319).sum(), (all_scores >= 11.319).mean())
262144 22 8.392333984375e-05
Now let us ask motif_pvalue() for the same
P-value:
um.motif_pvalue(ap1, score=11.319)
8.392333984375e-05
The two agree exactly. Listing every word quickly becomes impossible however,
since a motif of 20 positions already has about 1012 of them, and so
motif_pvalue() works differently. It builds the
full distribution of scores one column at a time, adding each column's scores
and their background probabilities to the distribution of the columns before
it. Because scores are kept on a grid of steps of 0.001, the distribution stays
small, and the result is exact for scores on that grid whatever the length of
the motif or the background. Several scores can be given at once:
um.motif_pvalue(ap1, score=[0, 5, 10, 13.356])
array([1.77230835e-02, 2.17819214e-03, 1.37329102e-04, 7.62939453e-06])
Note that even the maximum score has a P-value above 1/49, because two words reach it, one with a C at the centre and one with a G. Scores outside the range of the motif are fine too, and get P-values of one and zero, with a warning about the score that no word can reach:
with warnings.catch_warnings(record=True) as caught:
p = um.motif_pvalue(ap1, score=[-100, 100])
print(p)
print(caught[0].message)
[1. 0.] scores above the highest possible score get a P-value of 0 (motifs: AP-1)
From P-values to score thresholds¶
More often we want to go the other way, and find the score threshold for a
P-value cutoff. Given pvalue instead of score,
motif_pvalue() returns the lowest score on its
grid whose P-value is no greater than the cutoff:
cutoffs = [1e-2, 1e-3, 1e-4, 1e-5, 1e-6]
thresholds = um.motif_pvalue(ap1, pvalue=cutoffs)
thresholds
array([ 1.607, 7.126, 10.419, 13.04 , inf])
Let us check the threshold for 1e-4 by listing every word that reaches it with
get_matches(), and adding up their
probabilities:
matches = um.get_matches(ap1, thresholds[2])
print(len(matches), matches[:6])
um.prob_match(ap1, matches).sum()
26 ['TTGACTCAG', 'TTGAGTCAG', 'TTGACTCAA', 'TTGAGTCAA', 'TTGACTCAT', 'TTGAGTCAT']
np.float64(9.918212890625e-05)
These 26 words have a total probability just under 1e-4, as they should. The last two thresholds deserve a closer look. The threshold for 1e-5 is 13.04, which no word actually scores. Here are the four highest scores the motif can give:
np.unique(all_scores)[-4:]
array([12.448, 12.582, 13.039, 13.356])
Any threshold above 13.039 and up to 13.356 selects the same two words, those
with the maximum score, and 13.04 is the lowest such threshold on the grid. In
other words, a threshold need not be a score that any word can reach. The
threshold for 1e-6 is infinite, because even the best possible match, with its
P-value of about 7.6e-6, is not that unlikely. No word reaches an infinite
threshold, so get_matches() returns nothing and
scanning with such a cutoff finds no hits at all, which a warning points out
along with the smallest P-value the motif can give:
print(um.get_matches(ap1, thresholds[-1]))
targets = um.example_sequences()
with warnings.catch_warnings(record=True) as caught:
n = um.scan_sequences(ap1, targets, pvalue=1e-6).height
print(n)
print(caught[0].message)
[] 0 1 motif cannot reach a P-value of 1e-06 and will have no hits; the smallest P-value each can give is: AP-1 (7.63e-06)
This is worth keeping in mind with short or weak motifs. A cutoff as strict as 1e-6 is perfectly reasonable for a long motif, but most of the example motifs simply cannot reach it:
um.motif_pvalue(motifs, pvalue=1e-6)
array([13.312, inf, inf, inf, inf, inf, inf, inf,
15.395, 13.961, inf, inf])
With a list of motifs, motif_pvalue() gives
one value per motif. Only the TATA-box, ETS and CArG-box motifs, the three
10-position motifs, can produce a match with a P-value of 1e-6 or less.
Choosing a threshold¶
Let us now come back to thresholds given as fractions of the range of scores. For each example motif, here is the score at 80% of its range, the P-value of that score, and for comparison the threshold for a P-value of 1e-4:
rows = []
for m in motifs:
score80 = um.motif_score(m, 0.8)
rows.append({
"name": m.name,
"width": m.width,
"score_80pct": score80,
"pvalue_80pct": um.motif_pvalue(m, score=score80),
"threshold_1e-4": um.motif_pvalue(m, pvalue=1e-4),
})
pl.DataFrame(rows)
| name | width | score_80pct | pvalue_80pct | threshold_1e-4 |
|---|---|---|---|---|
| str | i64 | f64 | f64 | f64 |
| "TATA-box" | 10 | 4.649 | 0.003743 | 10.252 |
| "AP-1" | 9 | 3.6546 | 0.004848 | 10.419 |
| "CRE" | 8 | 3.788 | 0.003113 | 10.806 |
| "G-box" | 8 | 3.0718 | 0.009323 | 10.544 |
| "E-box" | 8 | 3.5714 | 0.00708 | 10.873 |
| "W-box" | 9 | 3.3868 | 0.004581 | 10.885 |
| "GATA" | 8 | 2.7588 | 0.014832 | 9.728 |
| "GC-box" | 9 | 5.2746 | 0.001308 | 10.158 |
| "ETS" | 10 | 4.2562 | 0.002893 | 10.596 |
| "CArG-box" | 10 | 3.0788 | 0.006936 | 10.544 |
| "Telo-box" | 9 | 4.7428 | 0.002068 | 10.277 |
| "CCAAT-box" | 9 | 4.117 | 0.006489 | 10.446 |
The same 80% threshold gives P-values ranging from about 0.0013 for the GC-box to about 0.015 for GATA, a more than ten-fold difference, which means that scanning with it would be far more lenient with some motifs than with others. A P-value cutoff on the other hand is equally strict for every motif, even though the scores it corresponds to vary.
Which P-value to choose is another matter. A useful way to think about it is the number of hits expected by chance: a 9-position motif has 992 windows in a sequence of 1000 bases, on each of two strands, so scanning the 100 example sequences at a cutoff of 1e-4 should turn up about 20 chance hits. The example sequences come as targets, with a W-box site planted in 60 of them, and controls with nothing planted. Let us scan both at a few cutoffs and count the sequences with at least one hit:
wbox = by_name["W-box"]
controls = um.example_sequences("controls")
def count_hits(motif, cutoffs):
rows = []
for p in cutoffs:
t = um.scan_sequences(motif, targets, pvalue=p)
c = um.scan_sequences(motif, controls, pvalue=p)
rows.append({
"pvalue": p,
"target_hits": t.height,
"target_seqs": t["sequence"].n_unique(),
"control_hits": c.height,
"control_seqs": c["sequence"].n_unique(),
})
return pl.DataFrame(rows)
count_hits(wbox, [1e-3, 1e-4, 1e-5])
| pvalue | 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, and the planted sites are lost among the chance matches. At 1e-4 the targets clearly stand out, though only 30 target sequences have a hit, so at least half of the 60 planted sites are missed; sites sampled from a motif are rarely perfect matches to it. At 1e-5 almost nothing is left in either set. So there is a trade-off between finding the real sites and keeping out the chance ones, and it is important to remember that there is likely never one single cutoff that works in every situation. It depends on the motif, on how much sequence is scanned, and on how many false positives can be tolerated.
Backgrounds¶
There is something odd about the controls in the table above: at 1e-4 they
have 10 hits, only half of the roughly 20 expected by chance. The reason is the
background. The background of a motif is used twice, once for the log ratios in
its PWM and once for the probability of each word, and the example motifs have
a uniform background while the example sequences are AT-rich, much like plant
promoters. get_bkg() shows by how much:
bkg = um.get_bkg(targets, k=1, rc=True)
at_rich = dict(zip(bkg["kmer"], bkg["probability"]))
bkg
| kmer | count | probability |
|---|---|---|
| str | i64 | f64 |
| "A" | 65804 | 0.32902 |
| "C" | 34196 | 0.17098 |
| "G" | 34196 | 0.17098 |
| "T" | 65804 | 0.32902 |
The W-box is fairly rich in G and C, which are uncommon in these sequences, so
its high-scoring words turn up less often than a uniform background predicts,
and its P-values are conservative for these sequences. The bkg_probs
argument of motif_pvalue() computes a P-value
or threshold as if the motif had a different background, which changes both
the PWM and the word probabilities. Setting the background of the motif itself
with replace() before it is turned into a PWM has the same effect:
wbox_at = wbox.replace(bkg=at_rich)
print(um.motif_pvalue(wbox, pvalue=1e-4))
print(um.motif_pvalue(wbox, pvalue=1e-4, bkg_probs=at_rich))
print(um.motif_pvalue(wbox_at, pvalue=1e-4))
10.885 10.599 10.599
Scanning uses the background stored in the motif, so let us scan again with the AT-rich version:
count_hits(wbox_at, [1e-4, 1e-5])
| pvalue | target_hits | target_seqs | control_hits | control_seqs |
|---|---|---|---|---|
| f64 | i64 | i64 | i64 | i64 |
| 0.0001 | 47 | 46 | 21 | 19 |
| 0.00001 | 12 | 12 | 0 | 0 |
Now the controls have 21 hits at 1e-4, very close to the 20 or so expected, so the P-values are much better calibrated for these sequences. At the same cutoffs more target sequences have a hit as well, 46 rather than 30 at 1e-4 and 12 rather than 1 at 1e-5. One last subtlety is that the background is only built into the scores when a motif is converted to a PWM. Changing the background of a motif that is already a PWM leaves its scores as they are, and changes only the probabilities of the words:
pwm = um.convert_type(ap1, "PWM")
pwm_at = pwm.replace(bkg=at_rich)
print(um.score_match(pwm, "ATGACTCAT"), um.score_match(pwm_at, "ATGACTCAT"))
print(um.prob_match(pwm, "ATGACTCAT"), um.prob_match(pwm_at, "ATGACTCAT"))
11.319 11.319 3.814697265625e-06 6.341174819652173e-06
So to move a motif to a new background, it is best to set the background on its counts or probabilities and convert to a PWM from there. The Choosing a background guide goes into how to pick one, with shuffled and composition-matched controls among the options.
Comparing columns¶
Let us now turn to comparing motifs. Motifs are compared one aligned column at
a time, and compare_motifs() uses the
Pearson correlation of the probabilities in each pair of columns. Two columns
can be compared directly with
compare_columns():
c1 = [0.7, 0.1, 0.1, 0.1]
c2 = [0.5, 0.0, 0.2, 0.3]
um.compare_columns(c1, c2, "PCC"), np.corrcoef(c1, c2)[0, 1]
(0.8006407690254359, np.float64(0.8006407690254357))
compare_columns() knows a number of other
metrics too, which can be handy for getting a feel for how they behave. Some
are distances, where lower means more alike, and others are similarities:
metrics = ["PCC", "WPCC", "EUCL", "WEUCL", "SEUCL", "MAN",
"HELL", "SW", "BHAT", "KL", "ALLR", "ALLR_LL"]
{m: round(um.compare_columns(c1, c2, m), 4) for m in metrics}
{'PCC': 0.8006,
'WPCC': 0.8006,
'EUCL': 0.3162,
'WEUCL': 0.2236,
'SEUCL': 0.1,
'MAN': 0.6,
'HELL': 0.3062,
'SW': 1.9,
'BHAT': 0.9062,
'KL': inf,
'ALLR': -inf,
'ALLR_LL': -2.0}
Note that KL and ALLR take logarithms of the probabilities, so the zero in the second column, facing a nonzero probability in the first, gives them an infinite score. Zeros in both columns at the same letter add nothing, but in general these two metrics need columns without zeros, which a small pseudocount takes care of:
c2_pseudo = [(x + 0.01) / 1.04 for x in c2]
um.compare_columns(c1, c2_pseudo, "KL"), um.compare_columns(c1, c2_pseudo, "ALLR")
(0.2871111380205665, 0.08999301862487857)
Comparing motifs¶
To compare two whole motifs, compare_motifs()
slides one along the other, and for DNA and RNA motifs along its reverse
complement too. At each offset it adds up the correlations of the overlapping
columns, and it keeps the alignment with the highest sum among those that
overlap by at least min_overlap columns. To see how this works, let us make
a few small motifs: one site, a shorter piece of it, the reverse complement of
a shifted version of it, and something unrelated.
site = um.create_motif("ATTGACCGAT", name="site", pseudocount=1, nsites=20)
part = um.create_motif("TTGACCG", name="part")
flipped = um.motif_rc(um.create_motif("GCTTGACCGA")).replace(name="flipped")
other = um.create_motif("CCCCCGGGCG", name="other")
small = [site, part, flipped, other]
print(flipped.consensus)
um.compare_motifs(small).round(3)
TCGGTCAAGC
array([[1. , 1. , 0.852, 0.333],
[1. , 1. , 1. , 0.556],
[0.852, 1. , 1. , 0.429],
[0.333, 0.556, 0.429, 1. ]])
By default the result is a matrix of the mean column correlation of the best
alignment of each pair, with the motifs as both rows and columns. part
matches site perfectly, whereas flipped only matches it once reverse
complemented, and even then its best score is 0.852 rather than 1. The reason
is that a shifted alignment overlaps by nine columns, eight of them identical
and one pairing an A with a C, and all overlapping columns count towards the
score. Without trying reverse complements, flipped looks rather unlike the
site, as the first row of the matrix shows:
um.compare_motifs(small, rc=False)[0].round(3)
array([1. , 1. , 0.259, 0.167])
To see the alignments themselves, give compare_to, one or more motifs to use
as queries, either as 0-based indices or as names. The result is then a
DataFrame with one row per pair, ordered by significance. Only pairs with a
q-value of at most qvalue are kept, so let us set it to one to see them all:
um.compare_motifs(small, compare_to="site", qvalue=1)
| subject | subject_i | target | target_i | offset | strand | overlap | score | subject_consensus | target_consensus | pvalue | qvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|
| str | i64 | str | i64 | i64 | str | i64 | f64 | str | str | f64 | f64 |
| "site" | 0 | "site" | 0 | 0 | "+" | 10 | 1.0 | "ATTGACCGAT" | "ATTGACCGAT" | 0.000007 | 0.000028 |
| "site" | 0 | "part" | 1 | 1 | "+" | 7 | 1.0 | "TTGACCG" | "TTGACCG" | 0.001013 | 0.001859 |
| "site" | 0 | "flipped" | 2 | 1 | "-" | 9 | 0.851852 | "ATTGACCGA" | "TCGGTCAAG" | 0.001394 | 0.001859 |
| "site" | 0 | "other" | 3 | -2 | "-" | 8 | 0.333333 | "TGACCGAT" | "CCCGGGCG" | 1.0 | 1.0 |
The offset is where the target starts relative to the query, strand is -
when the query was reverse complemented to line them up, and overlap is the
number of aligned columns. The two consensus columns show the aligned parts of
each motif, each in its own orientation, so on the - strand the reverse
complement of one lines up with the other: the reverse complement of TCGGTCAAG
is CTTGACCGA, which differs from ATTGACCGA only at its first position. The
score is the mean column correlation, the same as in the matrix.
min_overlap sets how short an alignment may be. Shorter overlaps let motifs
match by their edges, which can be what we want for motifs of quite different
lengths, but also makes it easier for unrelated motifs to find a few columns
that agree. For example, AP-1 and the CArG-box are aligned by only five columns
with the default of 5, so let us compare that with an overlap of at least 9:
carg = by_name["CArG-box"]
cols = ["target", "offset", "strand", "overlap", "score",
"subject_consensus", "target_consensus"]
for overlap in (5, 9):
hits = um.compare_motifs([ap1, carg], compare_to=0, qvalue=1, min_overlap=overlap)
print(hits.select(cols))
shape: (2, 7) āāāāāāāāāāāā¬āāāāāāāāā¬āāāāāāāāā¬āāāāāāāāāā¬āāāāāāāāāāā¬āāāāāāāāāāāāāāāāāāāā¬āāāāāāāāāāāāāāāāāāā ā target ā offset ā strand ā overlap ā score ā subject_consensus ā target_consensus ā ā --- ā --- ā --- ā --- ā --- ā --- ā --- ā ā str ā i64 ā str ā i64 ā f64 ā str ā str ā āāāāāāāāāāāāŖāāāāāāāāāŖāāāāāāāāāŖāāāāāāāāāāŖāāāāāāāāāāāŖāāāāāāāāāāāāāāāāāāāāŖāāāāāāāāāāāāāāāāāāā” ā AP-1 ā 0 ā + ā 9 ā 1.0 ā NTGASTCAN ā NTGASTCAN ā ā CArG-box ā -5 ā + ā 5 ā 0.335657 ā NTGAS ā WWTGG ā āāāāāāāāāāāā“āāāāāāāāā“āāāāāāāāā“āāāāāāāāāā“āāāāāāāāāāā“āāāāāāāāāāāāāāāāāāāā“āāāāāāāāāāāāāāāāāāā shape: (2, 7) āāāāāāāāāāāā¬āāāāāāāāā¬āāāāāāāāā¬āāāāāāāāāā¬āāāāāāāāāāāā¬āāāāāāāāāāāāāāāāāāāā¬āāāāāāāāāāāāāāāāāāā ā target ā offset ā strand ā overlap ā score ā subject_consensus ā target_consensus ā ā --- ā --- ā --- ā --- ā --- ā --- ā --- ā ā str ā i64 ā str ā i64 ā f64 ā str ā str ā āāāāāāāāāāāāŖāāāāāāāāāŖāāāāāāāāāŖāāāāāāāāāāŖāāāāāāāāāāāāŖāāāāāāāāāāāāāāāāāāāāŖāāāāāāāāāāāāāāāāāāā” ā AP-1 ā 0 ā + ā 9 ā 1.0 ā NTGASTCAN ā NTGASTCAN ā ā CArG-box ā 0 ā - ā 9 ā -0.026457 ā NTGASTCAN ā CCWWTWWTG ā āāāāāāāāāāāā“āāāāāāāāā“āāāāāāāāā“āāāāāāāāāā“āāāāāāāāāāāā“āāāāāāāāāāāāāāāāāāāā“āāāāāāāāāāāāāāāāāāā
Requiring nine columns forces the whole of AP-1 to be aligned, and the score
drops to around zero. A min_overlap longer than the shorter motif is the same
as asking for the shorter motif to be aligned in full.
Now let us compare all of the example motifs. The matrix is easier to read with the names alongside it:
names = [m.name for m in motifs]
scores = um.compare_motifs(motifs)
pl.DataFrame(scores.round(2), schema=names).insert_column(0, pl.Series("motif", names))
| motif | TATA-box | AP-1 | CRE | G-box | E-box | W-box | GATA | GC-box | ETS | CArG-box | Telo-box | CCAAT-box |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| "TATA-box" | 1.0 | 0.33 | 0.11 | 0.19 | 0.18 | 0.11 | 0.42 | -0.11 | 0.22 | 0.36 | 0.23 | 0.22 |
| "AP-1" | 0.33 | 1.0 | 0.46 | 0.25 | 0.23 | 0.81 | 0.25 | 0.3 | 0.22 | 0.34 | 0.28 | 0.38 |
| "CRE" | 0.11 | 0.46 | 1.0 | 0.58 | 0.29 | 0.5 | 0.18 | 0.17 | 0.45 | 0.2 | 0.19 | 0.39 |
| "G-box" | 0.19 | 0.25 | 0.58 | 1.0 | 0.52 | 0.46 | 0.37 | 0.32 | 0.51 | 0.23 | 0.19 | 0.22 |
| "E-box" | 0.18 | 0.23 | 0.29 | 0.52 | 1.0 | 0.53 | 0.49 | 0.21 | 0.32 | 0.2 | 0.21 | 0.3 |
| "W-box" | 0.11 | 0.81 | 0.5 | 0.46 | 0.53 | 1.0 | 0.18 | 0.33 | 0.25 | 0.2 | 0.53 | 0.44 |
| "GATA" | 0.42 | 0.25 | 0.18 | 0.37 | 0.49 | 0.18 | 1.0 | 0.18 | 0.46 | 0.41 | 0.48 | 0.58 |
| "GC-box" | -0.11 | 0.3 | 0.17 | 0.32 | 0.21 | 0.33 | 0.18 | 1.0 | 0.33 | -0.11 | 0.31 | 0.45 |
| "ETS" | 0.22 | 0.22 | 0.45 | 0.51 | 0.32 | 0.25 | 0.46 | 0.33 | 1.0 | 0.36 | 0.53 | 0.48 |
| "CArG-box" | 0.36 | 0.34 | 0.2 | 0.23 | 0.2 | 0.2 | 0.41 | -0.11 | 0.36 | 1.0 | 0.73 | 0.62 |
| "Telo-box" | 0.23 | 0.28 | 0.19 | 0.19 | 0.21 | 0.53 | 0.48 | 0.31 | 0.53 | 0.73 | 1.0 | 0.75 |
| "CCAAT-box" | 0.22 | 0.38 | 0.39 | 0.22 | 0.3 | 0.44 | 0.58 | 0.45 | 0.48 | 0.62 | 0.75 | 1.0 |
Note that the scores are symmetric, since swapping the query and target gives the same best alignment. The most similar pairs are AP-1 and the W-box at 0.81, followed by the Telo-box and the CCAAT-box at 0.75 and the Telo-box and the CArG-box at 0.73.
How significant is a comparison?¶
A mean correlation of 0.81 sounds high, but whether it is surprising depends on
how long the alignment is and on what columns in general look like. So
compare_motifs() also gives each alignment a
P-value: the probability of a sum of column correlations at least as high, for
an alignment of the same length, under a null model of random columns. This is
computed exactly from the distribution of column correlations, then
Bonferroni-corrected for the number of alignments that were tried. The
q-values are Benjamini-Hochberg adjusted P-values within each query. With the
default qvalue of 0.1, comparing AP-1 against everything else keeps just two
rows:
um.compare_motifs(motifs, compare_to="AP-1")
| subject | subject_i | target | target_i | offset | strand | overlap | score | subject_consensus | target_consensus | pvalue | qvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|
| str | i64 | str | i64 | i64 | str | i64 | f64 | str | str | f64 | f64 |
| "AP-1" | 1 | "AP-1" | 1 | 0 | "+" | 9 | 1.0 | "NTGASTCAN" | "NTGASTCAN" | 8.4899e-9 | 1.0188e-7 |
| "AP-1" | 1 | "W-box" | 5 | -2 | "+" | 7 | 0.808112 | "NTGASTC" | "TTGACYS" | 0.003279 | 0.019676 |
Apart from AP-1 itself, only the W-box passes. The two motifs share a TGAC,
which view_motifs() makes easy to see, since it
aligns the motifs the same way before drawing them:
fig = um.view_motifs([ap1, wbox])
The P-values and q-values are also available as matrices, with matrix_out
set to "pvalue" or "qvalue". Unlike the scores these need not be
symmetric, since both the null distribution and the adjustment are worked out
for each query. Let us look at the
q-values of the first six motifs:
qvalues = um.compare_motifs(motifs, matrix_out="qvalue")
first = names[:6]
pl.DataFrame(qvalues[:6, :6].round(4), schema=first).insert_column(0, pl.Series("motif", first))
| motif | TATA-box | AP-1 | CRE | G-box | E-box | W-box |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 |
| "TATA-box" | 0.0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 |
| "AP-1" | 1.0 | 0.0 | 1.0 | 1.0 | 1.0 | 0.0197 |
| "CRE" | 1.0 | 1.0 | 0.0001 | 0.3728 | 1.0 | 1.0 |
| "G-box" | 1.0 | 1.0 | 0.2687 | 0.0 | 0.5084 | 1.0 |
| "E-box" | 1.0 | 1.0 | 1.0 | 0.7721 | 0.0 | 1.0 |
| "W-box" | 1.0 | 0.0185 | 1.0 | 1.0 | 1.0 | 0.0 |
There are two null models to choose from. The default "empirical" null uses
the columns of all the motifs given, which describes what columns actually look
like but also means that the P-values depend on which motifs are compared.
With only AP-1 and the W-box, for example, the P-value of their alignment is
about twice as large:
um.compare_motifs([ap1, wbox], compare_to="AP-1", qvalue=1)
| subject | subject_i | target | target_i | offset | strand | overlap | score | subject_consensus | target_consensus | pvalue | qvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|
| str | i64 | str | i64 | i64 | str | i64 | f64 | str | str | f64 | f64 |
| "AP-1" | 0 | "AP-1" | 0 | 0 | "+" | 9 | 1.0 | "NTGASTCAN" | "NTGASTCAN" | 6.1252e-7 | 0.000001 |
| "AP-1" | 0 | "W-box" | 1 | -2 | "+" | 7 | 0.808112 | "NTGASTC" | "TTGACYS" | 0.006792 | 0.006792 |
The "parametric" null, which is for DNA and RNA motifs only, instead describes
random columns with a Dirichlet distribution around the background, given by
bkg and uniform by default. Its P-values do not depend on the other motifs,
which suits small sets of motifs where the empirical null has few columns to
go on:
cols = ["target", "score", "pvalue", "qvalue"]
print(um.compare_motifs(motifs, compare_to="AP-1", null="parametric").select(cols))
print(um.compare_motifs([ap1, wbox], compare_to="AP-1", null="parametric").select(cols))
shape: (2, 4) āāāāāāāāāā¬āāāāāāāāāāā¬āāāāāāāāāāāāā¬āāāāāāāāāāāāā ā target ā score ā pvalue ā qvalue ā ā --- ā --- ā --- ā --- ā ā str ā f64 ā f64 ā f64 ā āāāāāāāāāāŖāāāāāāāāāāāŖāāāāāāāāāāāāāŖāāāāāāāāāāāāā” ā AP-1 ā 1.0 ā 5.7428e-12 ā 6.8914e-11 ā ā W-box ā 0.808112 ā 0.001695 ā 0.010171 ā āāāāāāāāāā“āāāāāāāāāāā“āāāāāāāāāāāāā“āāāāāāāāāāāāā shape: (2, 4) āāāāāāāāāā¬āāāāāāāāāāā¬āāāāāāāāāāāāā¬āāāāāāāāāāāāā ā target ā score ā pvalue ā qvalue ā ā --- ā --- ā --- ā --- ā ā str ā f64 ā f64 ā f64 ā āāāāāāāāāāŖāāāāāāāāāāāŖāāāāāāāāāāāāāŖāāāāāāāāāāāāā” ā AP-1 ā 1.0 ā 5.7428e-12 ā 1.1486e-11 ā ā W-box ā 0.808112 ā 0.001695 ā 0.001695 ā āāāāāāāāāā“āāāāāāāāāāā“āāāāāāāāāāāāā“āāāāāāāāāāāāā
The parametric P-value for AP-1 against the W-box is the same in both sets, and the q-value changes only because there are fewer comparisons to adjust for. Note that a perfect match of a motif to itself is so unlikely under the parametric null that its P-value comes out as exactly zero.
Merging motifs¶
Motifs found in different experiments, or by different tools, are often variants of one another, and it can be useful to combine them into a single consensus motif. Let us make three variants of a W-box: two from slightly different sets of sites, one of them shifted by a position, and a third which is reverse complemented.
v1 = um.create_motif(
["CTTTGACCAA", "GTTTGACCAA", "CTTTGACTAA", "CTTTGACCAT"], name="W-box v1"
)
v2 = um.create_motif(["TTTGACCAAT", "TTTGACTAAT", "TTTGACCAAA"], name="W-box v2")
v3 = um.motif_rc(um.create_motif(
["ACTTTGACCA", "ACTTTGACCA", "AGTTTGACTA", "ACTTTGACCA", "TCTTTGACCA"]
)).replace(name="W-box v3")
variants = [v1, v2, v3]
[v.consensus for v in variants]
['CTTTGACCAA', 'TTTGACYAAW', 'TGGTCAAAGT']
merge_motifs() takes the motif with the highest
total information content as an anchor, aligns the others to it as
compare_motifs() would, reverse complementing
them where that fits better, and averages the aligned probabilities. The merged
motif spans every column that any of the motifs covers:
merged = um.merge_motifs(variants, new_name="W-box merged")
print(merged)
Motif name: W-box merged
Alphabet: DNA
Type: PPM
Strands: +-
Total IC: 20.16
Pseudocount: 0
Consensus: ACTTTGACYAAW
Target sites: 12
A C T T T G A C Y A A W
A 0.8 0.00 0 0 0 0 1 0 0.00 1 0.88 0.33
C 0.0 0.78 0 0 0 0 0 1 0.74 0 0.00 0.00
G 0.0 0.23 0 0 0 1 0 0 0.00 0 0.00 0.00
T 0.2 0.00 1 1 1 0 0 0 0.26 0 0.13 0.67
fig = um.view_motifs(variants + [merged])
By default every motif counts equally. With weighted=True the probabilities
are instead weighted by the number of sites behind each motif, which gives the
five sites of W-box v3 more say than the three of W-box v2. The differences
are small here, since the variants agree on most positions:
um.merge_motifs(variants, weighted=True).matrix.round(2)
array([[0.8 , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0. , 1. , 0.86,
0.33],
[0. , 0.78, 0. , 0. , 0. , 0. , 0. , 1. , 0.75, 0. , 0. ,
0. ],
[0. , 0.22, 0. , 0. , 0. , 1. , 0. , 0. , 0. , 0. , 0. ,
0. ],
[0.2 , 0. , 1. , 1. , 1. , 0. , 0. , 0. , 0.25, 0. , 0.14,
0.67]])
It is important to note that merge_motifs()
merges whatever it is given. It does not check whether the motifs are actually
alike, and an unrelated motif is merged in just the same, as long as it can
overlap the anchor by min_overlap columns. Adding the GC-box changes one
position of the consensus:
with_gc = um.merge_motifs(variants + [by_name["GC-box"]])
print(with_gc.name, with_gc.consensus)
merged_W-box v1+W-box v2+W-box v3+GC-box ACTTTGACCAAW
Deciding which motifs to merge is the job of
merge_similar(). It compares all the motifs,
links every pair with a q-value of at most qvalue in either direction, and
merges each connected group of motifs. Let us mix the variants in with the
example motifs, and first look at the clusters only:
mixed = motifs + variants
clusters = um.merge_similar(mixed, return_clusters=True)
clusters.select("name", "consensus", "motif_i", "cluster")
| name | consensus | motif_i | cluster |
|---|---|---|---|
| str | str | i64 | i64 |
| "TATA-box" | "KTATAWAWRN" | 0 | 0 |
| "AP-1" | "NTGASTCAN" | 1 | 1 |
| "CRE" | "TGACGTCA" | 2 | 2 |
| "G-box" | "NCACGTGN" | 3 | 3 |
| "E-box" | "NCAGCTGN" | 4 | 4 |
| "W-box" | "GBTTGACYS" | 5 | 5 |
| "GATA" | "CWGATART" | 6 | 6 |
| "GC-box" | "GGGGCGGGG" | 7 | 7 |
| "ETS" | "NCCGGAAGTG" | 8 | 8 |
| "CArG-box" | "CCWWTWWTGG" | 9 | 9 |
| "Telo-box" | "AAACCCTAA" | 10 | 10 |
| "CCAAT-box" | "CNCCAATNH" | 11 | 11 |
| "W-box v1" | "CTTTGACCAA" | 12 | 12 |
| "W-box v2" | "TTTGACYAAW" | 13 | 12 |
| "W-box v3" | "TGGTCAAAGT" | 14 | 12 |
The three variants form one cluster, and every other motif is on its own,
including the example W-box, whose sites are rather more varied. Without
return_clusters we get back one motif per cluster, with the variants
merged:
[m.name for m in um.merge_similar(mixed)]
['TATA-box', 'AP-1', 'CRE', 'G-box', 'E-box', 'W-box', 'GATA', 'GC-box', 'ETS', 'CArG-box', 'Telo-box', 'CCAAT-box', 'merged_W-box v1+W-box v2+W-box v3']
The default cutoff of 0.01 is fairly strict. Loosening it pulls more motifs together, and at 0.1 one cluster gathers six of them:
loose = um.merge_similar(mixed, qvalue=0.1, return_clusters=True)
loose.filter(pl.col("cluster") == 1).select("name", "consensus", "cluster")
| name | consensus | cluster |
|---|---|---|
| str | str | i64 |
| "AP-1" | "NTGASTCAN" | 1 |
| "W-box" | "GBTTGACYS" | 1 |
| "CArG-box" | "CCWWTWWTGG" | 1 |
| "W-box v1" | "CTTTGACCAA" | 1 |
| "W-box v2" | "TTTGACYAAW" | 1 |
| "W-box v3" | "TGGTCAAAGT" | 1 |
To see why, we can list the pairs that pass the same cutoff with
compare_motifs(), leaving out the
comparisons of each motif with itself:
links = um.compare_motifs(mixed, compare_to=list(range(len(mixed))), qvalue=0.1)
links.filter(pl.col("subject") != pl.col("target")).select(
"subject", "target", "overlap", "score", "qvalue"
)
| subject | target | overlap | score | qvalue |
|---|---|---|---|---|
| str | str | i64 | f64 | f64 |
| "W-box v1" | "W-box v3" | 9 | 0.999137 | 1.1796e-7 |
| "W-box v3" | "W-box v1" | 9 | 0.999137 | 7.4148e-7 |
| "W-box v1" | "W-box v2" | 9 | 0.991949 | 0.000014 |
| "W-box v2" | "W-box v1" | 9 | 0.991949 | 0.000016 |
| "W-box v2" | "W-box v3" | 8 | 0.995698 | 0.000026 |
| "W-box v3" | "W-box v2" | 8 | 0.995698 | 0.000129 |
| "W-box" | "AP-1" | 7 | 0.808112 | 0.025175 |
| "AP-1" | "W-box" | 7 | 0.808112 | 0.027492 |
| "W-box" | "W-box v2" | 8 | 0.689935 | 0.05212 |
| "W-box" | "W-box v1" | 9 | 0.62429 | 0.065875 |
| "W-box" | "W-box v3" | 9 | 0.612002 | 0.067118 |
| "W-box v2" | "W-box" | 8 | 0.689935 | 0.067483 |
| "W-box v2" | "CArG-box" | 5 | 0.88055 | 0.069788 |
The variants are linked to the example W-box, which is linked to AP-1 through
their shared TGAC, and W-box v2 is also linked to the CArG-box by a
five-column overlap of their AT-rich ends. Since clusters are connected groups,
a single link is enough to join two motifs to the same cluster, however unlike
they are otherwise. Whether that is desirable depends on what the motifs are
for, so it is always worth looking at the clusters before trusting the merged
motifs.
Motif trees¶
A tree is often the clearest way of seeing how a set of motifs relate to one
another. motif_tree() compares the motifs with
compare_motifs(), turns the mean column
correlations into distances as (1 - score) / 2, so that identical motifs are
at a distance of zero and perfectly anticorrelated ones at one, and clusters
the distances hierarchically. The default is complete linkage, and method
offers the others:
tree = um.motif_tree(motifs)
tree
MotifTree(labels=['TATA-box', 'AP-1', 'CRE', 'G-box', 'E-box', 'W-box', 'GATA', 'GC-box', 'ETS', 'CArG-box', 'Telo-box', 'CCAAT-box'], method='complete')
Note that min_overlap is 6 here by default, rather than the 5 of
compare_motifs(). The returned
MotifTree holds the clustering in a few forms. The
tip labels come in input order, order gives the motifs in the order they
appear down the tree, and height the height of each merge:
print([tree.labels[i] for i in tree.order])
tree.height.round(3)
['E-box', 'GATA', 'TATA-box', 'CArG-box', 'Telo-box', 'CCAAT-box', 'ETS', 'CRE', 'G-box', 'GC-box', 'AP-1', 'W-box']
array([0.096, 0.125, 0.204, 0.211, 0.254, 0.273, 0.35 , 0.389, 0.411,
0.417, 0.579])
The tree is also available as a Newick string, ready for any tree viewer, and
MotifTree.write_newick() saves it to
a file:
print(tree.newick)
import tempfile, pathlib
newick_path = pathlib.Path(tempfile.mkdtemp()) / "motifs.nwk"
tree.write_newick(newick_path)
(((E-box:0.1270868944,GATA:0.1270868944):0.0782601072,(TATA-box:0.1946082618,(CArG-box:0.1017729637,(Telo-box:0.06233023255,CCAAT-box:0.06233023255):0.03944273116):0.09283529808):0.01073873983):0.08398608608,((ETS:0.1366146905,(CRE:0.1053002703,G-box:0.1053002703):0.03131442018):0.07178003125,(GC-box:0.1749309699,(AP-1:0.04797203635,W-box:0.04797203635):0.1269589336):0.03346375184):0.08093836595);
For further analysis, linkage gives the same tree as a SciPy linkage matrix,
so all of the scipy.cluster.hierarchy functions can be used with it, such
as fcluster() to cut it into a set number of clusters:
from scipy.cluster import hierarchy
hierarchy.fcluster(tree.linkage, 4, criterion="maxclust")
array([2, 4, 3, 3, 1, 4, 1, 4, 3, 2, 2, 2], dtype=int32)
Cutting a tree is also built in, as
MotifTree.cutree(), which takes either a
number of groups k or a height h to cut at. Groups are numbered from zero
in the order they first appear among the motifs:
groups = tree.cutree(k=4)
pl.DataFrame({"motif": names, "group": groups})
| motif | group |
|---|---|
| str | i64 |
| "TATA-box" | 0 |
| "AP-1" | 1 |
| "CRE" | 2 |
| "G-box" | 2 |
| "E-box" | 3 |
| "W-box" | 1 |
| "GATA" | 3 |
| "GC-box" | 1 |
| "ETS" | 2 |
| "CArG-box" | 0 |
| "Telo-box" | 0 |
| "CCAAT-box" | 0 |
The four groups are the same as those from fcluster(), only numbered
differently. Cutting lower down, at a height of 0.3, splits them further:
tree.cutree(h=0.3)
array([0, 1, 2, 2, 3, 1, 3, 4, 2, 5, 5, 5])
MotifTree.plot() draws the tree with
matplotlib. With tip_logos=True each motif's logo is drawn at its tip,
aligned to the others as in view_motifs(), which
makes it much easier to see what drives the clustering. Branches can be
coloured by group, given either as one value per motif or as a motif
attribute, and a branch takes the colour of a group when every motif below it
belongs to that group. Let us colour the branches by the four groups
from MotifTree.cutree(), leaving out the
legend since the groups have no names:
ax = tree.plot(groups=groups, tip_logos=True, legend=False)
Tip logos need the rectangular or slanted layouts. The circular and fan
layouts suit larger trees, and tip_size can size the tips by a motif
attribute such as the total information content. Let us colour the branches by
family this time. The colours are taken in turn from the matplotlib colour
cycle, in the order the families first appear, so the TBP family of the
TATA-box is blue and the bZIP family orange; colours sets them explicitly
instead. The legends beside the tree name the families and show the tip sizes.
ax = tree.plot(layout="circular", groups="family", tip_size="icscore")
The three bZIP motifs do not all end up together. CRE and the G-box, which share an ACGT core, sit next to each other, but AP-1 is closest to the W-box, as we saw earlier.
Finally, motif_tree() also accepts a square
matrix of distances instead of motifs, for trees built from some other measure
of similarity. In that case the labels have to be given, and since there are
no motifs, there can be no tip logos. Building the same distances as before by
hand gives the same tree:
distances = (1 - um.compare_motifs(motifs, min_overlap=6)) / 2
same_tree = um.motif_tree(distances, labels=names)
same_tree.newick == tree.newick
True
Where next¶
The API reference describes every argument of the functions used here, in its P-values and scores and comparison and clustering pages. The P-values of individual matches are what scanning and enrichment testing are built on, and the scanning and enrichment pages cover those.