Sequence motifs¶
Motifs are a more practical representation of the short, recurring elements found in biological sequences than a single consensus string, since they describe how variable each position of a site is. There are a few ways of representing a motif as a matrix: as counts, as probabilities, as log-odds weights, as information content, and as contribution weights from attribution analyses. Naming conventions vary somewhat between tools, but in this package these are the PCM, PPM, PWM, ICM and CWM types.
This guide goes through each of them, how they are calculated from one another
and what they are good for, using a single small motif throughout. Turning one
type into another is always a call to
convert_type(), and the
Working with motifs guide covers the rest of what
can be done with a Motif.
import numpy as np
import universalmotif as um
Position count matrices¶
Also known as position frequency matrices, these are the most basic representation of a motif. Simply, for each position the number of times each letter was seen is counted. For example, let us start from the following six sites:
sites = ["AAGAAT", "ATCATA", "AAGTAA", "AACAAA", "ATTAAA", "AAGAAT"]
pcm = um.create_motif(sites, type="PCM", name="example")
print(pcm)
Motif name: example
Alphabet: DNA
Type: PCM
Strands: +-
Total IC: 7.4
Pseudocount: 0
Consensus: AWSAAW
Target sites: 6
A W S A A W
A 6 4 0 5 5 4
C 0 0 2 0 0 0
G 0 0 3 0 0 0
T 0 2 1 1 1 2
Each column of the matrix is a position and each row is a letter, and the
column names are the consensus letters, which are covered further down. Note
that every position adds up to six, the number of sites the motif was made
from, which is kept as nsites. The matrix itself is a numpy array with one
row per letter.
print(pcm.letters, pcm.nsites)
pcm.matrix.sum(axis=0)
['A', 'C', 'G', 'T'] 6.0
array([6., 6., 6., 6., 6., 6.])
Position probability matrices¶
Also sometimes known as position frequency matrices or position weight matrices, these hold the probability of each letter at each position. The probability P of a letter N at a single position is calculated from the counts C as:
In numpy this is a matter of dividing each column by its sum, which gives the same matrix as converting the motif to the PPM type:
counts = pcm.matrix
print((counts / counts.sum(axis=0)).round(3))
ppm = um.convert_type(pcm, "PPM")
ppm.matrix.round(3)
[[1. 0.667 0. 0.833 0.833 0.667] [0. 0. 0.333 0. 0. 0. ] [0. 0. 0.5 0. 0. 0. ] [0. 0.333 0.167 0.167 0.167 0.333]]
array([[1. , 0.667, 0. , 0.833, 0.833, 0.667],
[0. , 0. , 0.333, 0. , 0. , 0. ],
[0. , 0. , 0.5 , 0. , 0. , 0. ],
[0. , 0.333, 0.167, 0.167, 0.167, 0.333]])
Every position now sums to 1. From this representation, the probability of any combination of letters can be calculated by multiplying the probabilities of its letters at each position. For example, AAGAAA has a probability of about 15%:
def site_probability(motif, site):
rows = [motif.letters.index(letter) for letter in site]
return np.prod(motif.matrix[rows, range(len(site))])
site_probability(ppm, "AAGAAA"), site_probability(ppm, "AAAAAA")
(np.float64(0.154320987654321), np.float64(0.0))
However, when a motif is made from only a few sites, many zeroes appear in the matrix. This means that the probability of AAAAAA is currently zero, since no site had an A at the third position. When scanning large amounts of sequence, ruling out combinations such as these is usually undesirable, as mismatches do occur in real binding sites.
Pseudocounts¶
This is fixed by adding a pseudocount p to the counts. Usually a small number such as 1, it is introduced into the calculation as such:
The pseudocount is added to the top and bottom of the fraction, but on top,
where it is specific to each letter, it is divided by the number of letters
n in the alphabet, which is 4 for DNA. Every motif carries its own
pseudocount which is used whenever it is converted from one type to another,
and convert_type() can also be given one
directly:
ppm1 = um.convert_type(pcm, "PPM", pseudocount=1)
ppm1.matrix.round(3)
array([[0.893, 0.607, 0.036, 0.75 , 0.75 , 0.607],
[0.036, 0.036, 0.321, 0.036, 0.036, 0.036],
[0.036, 0.036, 0.464, 0.036, 0.036, 0.036],
[0.036, 0.321, 0.179, 0.179, 0.179, 0.321]])
Now, though unlikely, AAAAAA is no longer impossible.
site_probability(ppm1, "AAGAAA"), site_probability(ppm1, "AAAAAA")
(np.float64(0.08595509114041022), np.float64(0.006611930087723863))
Since there are only six sites here, a pseudocount of 1 has a large effect on the probabilities. This effect shrinks as the number of sites grows, which we can see by building the same motif from sixty sites instead, with each of the six sites repeated ten times. The probabilities without a pseudocount are identical, yet the pseudocount barely moves them.
pcm60 = um.create_motif(sites * 10, type="PCM")
um.convert_type(pcm60, "PPM", pseudocount=1).matrix.round(3)
array([[0.988, 0.66 , 0.004, 0.824, 0.824, 0.66 ],
[0.004, 0.004, 0.332, 0.004, 0.004, 0.004],
[0.004, 0.004, 0.496, 0.004, 0.004, 0.004],
[0.004, 0.332, 0.168, 0.168, 0.168, 0.332]])
The pseudocount is set when a motif is created, or changed later with
replace(). Note that changing it does not touch the matrix of a PPM; it is
only applied the next time the motif is converted to another type, and
normalize() can be used to apply it straight
away.
Position weight matrices¶
The position weight matrix, also known as the position-specific scoring matrix or log-odds matrix, gives every letter at every position a score which can be used to evaluate how well a sequence matches the motif. Though there are several ways of calculating these scores S, the most common is the log of each probability relative to the background frequency B of that letter:
Letters that are more likely in the motif than in the background get positive scores, and letters that are less likely get negative ones. Starting from the PPM without a pseudocount and a uniform background, where each letter has a probability of 0.25, this gives:
pwm0 = um.convert_type(pcm, "PWM")
pwm0.matrix.round(3)
array([[ 2. , 1.415, -inf, 1.737, 1.737, 1.415],
[ -inf, -inf, 0.415, -inf, -inf, -inf],
[ -inf, -inf, 1. , -inf, -inf, -inf],
[ -inf, 0.415, -0.585, -0.585, -0.585, 0.415]])
Note that the positions no longer have equal sums. To score a sequence, the
scores for its letters at each position are added up, which
score_match() does for us. AAGAAA has a score
of 9.302, but much like with the PPM, AAAAAA can never be recovered with this
motif since its score is minus infinity. By default, a motif whose scores
would be infinite is given a pseudocount of 1 before scoring, so we need to ask
for infinite scores explicitly:
um.score_match(pwm0, ["AAGAAA", "AAAAAA"], allow_nonfinite=True)
array([9.302, -inf])
This can be avoided simply by starting from the pseudocount-adjusted probabilities, so let us give the motif a pseudocount of 1 and convert it again:
pcm = pcm.replace(pseudocount=1)
pwm = um.convert_type(pcm, "PWM")
pwm.matrix.round(3)
array([[ 1.837, 1.28 , -2.807, 1.585, 1.585, 1.28 ],
[-2.807, -2.807, 0.363, -2.807, -2.807, -2.807],
[-2.807, -2.807, 0.893, -2.807, -2.807, -2.807],
[-2.807, 0.363, -0.485, -0.485, -0.485, 0.363]])
um.score_match(pwm, ["AAGAAA", "AAAAAA"])
array([8.457, 4.757])
Now, the score for AAGAAA is 8.456, and the score for AAAAAA is 4.757. Though the latter is low, it is no longer minus infinity because of a single mismatch. (Note that scores are truncated to three decimal places, which keeps them consistent with the P-values calculated by the package.)
When searching for this motif, one more consideration is needed, which is a
minimum score. For example, AAGAAA and ATTTTT both have positive scores, but
one is much higher than the other. Should ATTTTT then be discarded? To answer
this, a threshold is set, often as a fraction of the highest possible score.
motif_score() gives the lowest and highest
possible scores, or the score at any fraction of the range:
print(um.motif_score(pwm))
print(um.motif_score(pwm, 0.25, threshold_type="fromzero"))
um.score_match(pwm, "ATTTTT")
[-16.842 8.457] 2.11425
1.105
The highest possible score is 8.457, so a threshold of 25% of it gives a minimum score of 2.114, and ATTTTT with its score of 1.105 would indeed be discarded.
Another way of choosing a minimum score is to start from a P-value instead,
which is the probability of a random sequence scoring at least as high.
motif_pvalue() calculates these exactly in
both directions, and they are what
scan_sequences() uses for its cutoffs:
print(um.motif_pvalue(pwm, score=4.757))
um.motif_pvalue(pwm, pvalue=0.001)
0.006591796875
7.08
So a score of 4.757, which is what AAAAAA gets, is reached or beaten by about 0.66% of random sites, while a P-value of 0.001 needs a score of at least 7.08. The P-values and scores reference has more on these functions.
Whichever method you use however, it is important to remember that there is likely never one single threshold that works in every situation. Sequences passing a lenient threshold can still be bound by transcription factors that bind DNA with less specificity or are present at higher concentrations, and the opposite can be true of more specific or less abundant ones.
Backgrounds¶
So far the background has been uniform, but real sequences rarely are. Every
motif has a background, bkg, which is used for PWMs and, as we will see
further down, for one way of calculating information content. Let us suppose
our motif comes from AT-rich sequences, where A and T each make up 40% of the
letters:
at_rich = {"A": 0.4, "C": 0.1, "G": 0.1, "T": 0.4}
pcm_at = pcm.replace(bkg=at_rich)
um.convert_type(pcm_at, "PWM").matrix.round(3)
array([[ 1.158, 0.602, -3.485, 0.907, 0.907, 0.602],
[-1.485, -1.485, 1.684, -1.485, -1.485, -1.485],
[-1.485, -1.485, 2.215, -1.485, -1.485, -1.485],
[-3.485, -0.316, -1.163, -1.163, -1.163, -0.316]])
The scores for A and T have dropped, since these letters are no longer as
surprising, and the G at position 3 has gone from 0.893 to 2.215 now that G is
rare in the background. The backgrounds of real sequences can be measured
with get_bkg().
Information content matrices¶
Finally, the information content matrix. This type aims to include another
consideration, which is whether some positions are more important than others.
To explore this, let us draw the PPM as a sequence logo with
view_motifs():
fig = um.view_motifs(ppm, use_type="PPM")
In this case, the height of each letter is its probability, so that every position has the same height. However, when drawn as an ICM, the default for logos, the motif looks like this:
fig = um.view_motifs(ppm)
Now, the total height of each position is scaled by the information at that position. Simply put, the information at a position is an indication of how conserved it is. In this motif the first position is highly conserved, always being an A, whereas the third position is much less so, as no single letter has a high probability there.
For every position, the total height is calculated by taking the maximum possible information content, or maximum uncertainty, and subtracting the actual uncertainty at that position. This is based on Shannon's entropy, with the final values in bits. The maximum depends on the number of letters n in the alphabet:
For DNA this is 2 bits. The actual uncertainty U of a position is calculated as:
Here, the uncertainty of the position is the sum of the uncertainties contributed by each letter. The information content of the position is then:
and the height of each letter is its probability multiplied by this:
Let us do this calculation with numpy and compare it to converting the motif
to the ICM type, keeping in mind that 0 * log2(0) counts as zero:
P = ppm.matrix
with np.errstate(divide="ignore", invalid="ignore"):
U = -np.nansum(P * np.log2(P), axis=0)
IC = P * (np.log2(4) - U)
print(IC.round(3))
um.convert_type(ppm, "ICM").matrix.round(3)
[[2. 0.721 0. 1.125 1.125 0.721] [0. 0. 0.18 0. 0. 0. ] [0. 0. 0.27 0. 0. 0. ] [0. 0.361 0.09 0.225 0.225 0.361]]
array([[2. , 0.721, 0. , 1.125, 1.125, 0.721],
[0. , 0. , 0.18 , 0. , 0. , 0. ],
[0. , 0. , 0.27 , 0. , 0. , 0. ],
[0. , 0.361, 0.09 , 0.225, 0.225, 0.361]])
Note that no position can add up to more than 2. The information content of a
single position is also available from
position_icscore(), here calculated from
the counts without a pseudocount, and the sum over all positions is the
motif's icscore, which is the total IC shown when printing a motif.
per_position = [um.position_icscore(col, type="PCM", pseudocount=0) for col in pcm.matrix.T]
print(np.round(per_position, 3))
print(round(sum(per_position), 3), round(ppm.icscore, 3))
[2. 1.082 0.541 1.35 1.35 1.082] 7.404 7.404
In its original formulation this calculation also includes a correction for small numbers of sites, which is rarely used and is not applied here.
Relative entropy¶
An alternative to Shannon's entropy is relative entropy, also known as the Kullback-Leibler divergence. While the information content above has the advantage of the same maximum at every position, it takes no account of the background. Relative entropy does, at the cost of positions no longer sharing the same maximum:
This can give values below zero for letters rarer in the motif than in the
background, which are replaced with zero. Pass relative_entropy=True to
convert_type() to use it. With a uniform
background the motif becomes:
um.convert_type(ppm, "ICM", relative_entropy=True).matrix.round(3)
array([[2. , 0.943, 0. , 1.447, 1.447, 0.943],
[0. , 0. , 0.138, 0. , 0. , 0. ],
[0. , 0. , 0.5 , 0. , 0. , 0. ],
[0. , 0.138, 0. , 0. , 0. , 0.138]])
The motif looks rather different with the AT-rich background from earlier:
ppm_at = ppm.replace(bkg=at_rich)
um.convert_type(ppm_at, "ICM", relative_entropy=True).matrix.round(3)
array([[1.322, 0.491, 0. , 0.882, 0.882, 0.491],
[0. , 0. , 0.579, 0. , 0. , 0. ],
[0. , 0. , 1.161, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0. ]])
Notice how the height of A at position 1 has dropped from 2 to 1.322, and how the height of G at position 3 has risen from 0.5 to 1.161. This makes sense when you consider that in these sequences you are more likely to get an A and less likely to get a G to begin with. Logos can be drawn the same way, and here we use the counts with the AT-rich background from earlier, which have a pseudocount of 1.
fig = um.view_motifs(pcm_at, relative_entropy=True)
Consensus letters¶
The column names of a printed motif are its consensus letters, which use the IUPAC ambiguity codes for DNA and RNA so that positions without a single dominant letter can still be described. For example, W stands for A or T, S for C or G, and N for any letter. Our motif's consensus is AWSAAW:
ppm.consensus
'AWSAAW'
The letter for a single position comes from
get_consensus(). A letter is called on its
own when it clearly dominates the position, a pair or trio of letters when
they share most of the probability between them, and N otherwise:
for column in ([0.7, 0.1, 0.1, 0.1], [0.45, 0.05, 0.05, 0.45],
[0.3, 0.3, 0.3, 0.1], [0.25, 0.25, 0.25, 0.25]):
print(column, um.get_consensus(column))
[0.7, 0.1, 0.1, 0.1] A [0.45, 0.05, 0.05, 0.45] W [0.3, 0.3, 0.3, 0.1] V [0.25, 0.25, 0.25, 0.25] N
Going the other way, consensus_to_ppm()
gives the probabilities for a consensus letter, keeping a tiny probability for
the letters it excludes. create_motif() shares
out the probabilities in the same way when making a motif from a consensus
string, though with exact zeroes, which is why S is split evenly between C and
G here:
print(um.consensus_to_ppm("S"))
print(um.create_motif("CCNSNGG"))
[0.001 0.499 0.499 0.001]
Motif name: motif
Alphabet: DNA
Type: PPM
Strands: +-
Total IC: 9
Pseudocount: 0
Consensus: CCNSNGG
C C N S N G G
A 0 0 0.25 0.0 0.25 0 0
C 1 1 0.25 0.5 0.25 0 0
G 0 0 0.25 0.5 0.25 1 1
T 0 0 0.25 0.0 0.25 0 0
Amino acid motifs have their own consensus letters, which come from
get_consensus_aa().
Contribution weight matrices¶
A contribution weight matrix (CWM) stores signed, unnormalized contributions from an attribution analysis of a trained sequence model. Positive and negative entries represent effects in opposite directions under that analysis, and they are neither probabilities nor log-odds scores. How to interpret them depends on the model, the prediction task, the attribution method and the reference used to calculate them.
Since a signed matrix otherwise looks like a PWM, a CWM must always be asked
for with type="CWM":
contributions = np.array([[1.0, -0.25],
[-0.5, 1.0],
[0.0, -0.5],
[-0.5, 0.0]])
cwm = um.create_motif(contributions, alphabet="DNA", type="CWM", name="contributions")
fig = um.view_motifs(cwm, use_type="CWM")
Letters with negative contributions are drawn below the line. CWMs can be read
from files with cwm=True in read_meme() or
type="CWM" in read_matrix(), and
trim_cwm() trims weak positions from their edges
using absolute contributions. Note that the number of sites is unknown unless
given, since a column sum is not a number of contributing sequences.
Converting a CWM to a PPM normalizes the absolute values of each column, which throws away both their signs and their scale.
um.convert_type(cwm, "PPM").matrix.round(3)
array([[0.5 , 0.143],
[0.25 , 0.571],
[0. , 0.286],
[0.25 , 0. ]])
This gives a probability matrix which other functions can work with, but it does not recover a binding model, and P-values from scanning with it are not probabilities of anything the model learnt. There is also no way back from a PPM to the original CWM, so it is worth keeping the CWM alongside anything derived from it.
Logos for each type¶
Since the types hold such different values, their logos look quite different too. We have already seen the PPM, ICM and CWM; here is the same motif as counts and as log-odds weights, where letters with negative scores are drawn below the line.
fig = um.view_motifs(pcm, use_type="PCM")
fig = um.view_motifs(pcm, use_type="PWM")
Counts logos only differ from probability logos in their scale, and weight logos tend to be dominated by the large negative scores of letters that were never seen. For most purposes the ICM logo is the most informative, which is why it is the default. Note that drawing a PWM needs finite scores, so a motif without a pseudocount cannot be drawn this way.
Where next¶
The Working with motifs guide shows how to create, read, write and modify motifs, and the API reference describes everything else.