Gapped motifs and other alphabets¶
Most motifs are short, contiguous stretches of DNA, and most of the package is written with those in mind. Not every motif fits that mould, however. Some binding sites are made of two conserved halves separated by a spacer of variable length, and some motifs are not made of DNA at all, but of RNA, amino acids, or letters of your own choosing. This guide shows how to handle both: first gapped motifs, then motifs in other alphabets.
As usual, everything comes from the top-level package:
import warnings
import polars as pl
import universalmotif as um
Adding gaps to a motif¶
Let us consider a motif made of two copies of the half-site AGGTCA. Proteins
which bind such direct repeats often tolerate a spacer of a few bases between
the two halves, and how long that spacer is can vary from one site to the
next. An ordinary motif can only describe one spacing at a time, so to
describe all of them at once we can add a gap with
add_gap(). The gap comes after the first
gaploc positions, and can span anywhere from mingap to maxgap
positions. Here the gap sits between the two half-sites and spans one to five
positions:
dr = um.create_motif("AGGTCAAGGTCA", name="DR", nsites=20, pseudocount=1)
dr_gapped = um.add_gap(dr, gaploc=6, mingap=1, maxgap=5)
print(dr_gapped)
Motif name: DR
Alphabet: DNA
Type: PPM
Strands: +-
Total IC: 20.65
Pseudocount: 1
Consensus: AGGTCA..AGGTCA
Target sites: 20
Gap locations: 6-7
Gap sizes: 1-5
A G G T C A A G G T C A
A 0.96 0.01 0.01 0.01 0.01 0.96 0.96 0.01 0.01 0.01 0.01 0.96
C 0.01 0.01 0.01 0.01 0.96 0.01 0.01 0.01 0.01 0.01 0.96 0.01
G 0.01 0.96 0.96 0.01 0.01 0.01 0.01 0.96 0.96 0.01 0.01 0.01
T 0.01 0.01 0.01 0.96 0.01 0.01 0.01 0.01 0.01 0.96 0.01 0.01
Note that the matrix itself is unchanged, since the gap has no columns of its own. The printed consensus marks where the gap goes with two dots whatever its size, and the gap's location and range of sizes are listed above the matrix. The same information is available as attributes.
print(dr_gapped.isgapped, dr_gapped.gaploc, dr_gapped.mingap, dr_gapped.maxgap)
True [6] [1] [5]
A motif can have more than one gap, in which case gaploc, mingap and
maxgap take one value per gap, with shorter lists being recycled. A gap may
also be allowed to have a size of zero, which means the motif can be matched
with or without it. For example, here is a version of the same motif with an
additional optional gap after the third position:
dr_two_gaps = um.add_gap(dr, gaploc=[3, 6], mingap=[0, 1], maxgap=[1, 5])
print(dr_two_gaps.gaploc, dr_two_gaps.mingap, dr_two_gaps.maxgap)
[3, 6] [0, 1] [1, 5]
This motif can be matched in 2 × 5 = 10 different ways, one for each combination of gap sizes, which will become important later when we look at P-values.
Scanning gapped motifs¶
Now, let us make a few random sequences and plant a direct repeat in each, with a spacer of one, two and four bases in the first three. The fourth gets a repeat with a three base spacer, planted as its reverse complement so that it can only be found on the reverse strand:
def plant(seq, site, at=90):
return seq[:at] + site + seq[at + len(site):]
sites = {
"seq1": "AGGTCA" + "A" + "AGGTCA",
"seq2": "AGGTCA" + "GA" + "AGGTCA",
"seq3": "AGGTCA" + "CATG" + "AGGTCA",
"seq4": "TGACCT" + "GAC" + "TGACCT",
}
random_seqs = um.create_sequences(seqnum=4, seqlen=200, seed=4)
seqs = {name: plant(s, site) for (name, site), s in zip(sites.items(), random_seqs)}
The original motif, which has no gap, expects the two half-sites to sit right next to each other, and so it finds none of them.
um.scan_sequences(dr, seqs, pvalue=1e-5).height
0
The gapped motif, on the other hand, finds all four. When given a gapped
motif, scan_sequences() tries every
combination of gap sizes at every position, and shows the letters skipped
over by the gap as dots in the match column:
cols = ["sequence", "start", "end", "strand", "score", "match", "pvalue"]
hits = um.scan_sequences(dr_gapped, seqs, pvalue=1e-5)
hits.select(cols)
| sequence | start | end | strand | score | match | pvalue |
|---|---|---|---|---|---|---|
| str | i64 | i64 | str | f64 | str | f64 |
| "seq1" | 90 | 103 | "+" | 23.364 | "AGGTCA.AGGTCA" | 5.9605e-8 |
| "seq2" | 90 | 104 | "+" | 23.364 | "AGGTCA..AGGTCA" | 5.9605e-8 |
| "seq3" | 90 | 106 | "+" | 23.364 | "AGGTCA....AGGTCA" | 5.9605e-8 |
| "seq4" | 90 | 105 | "-" | 23.364 | "AGGTCA...AGGTCA" | 5.9605e-8 |
Gapped motifs are scanned on both strands like any other DNA motif, so the
site in seq4 is found on the - strand, with its match shown in the
orientation of the motif. Note also that start and end include the gap,
which is why the hits are of different lengths.
The P-value of a gapped hit is simply that of its score under the motif
without its gaps. Since the letters inside the gap are not scored at all,
every hit here has the same score and P-value, and we can get the very same
P-value from motif_pvalue() using the
ungapped motif.
um.motif_pvalue(dr, score=hits["score"][0])
5.960464477539063e-08
If you would rather scan a gapped motif as if it had no gaps, set
use_gaps=False. Here, that gives the same result as scanning the original
motif.
um.scan_sequences(dr_gapped, seqs, pvalue=1e-5, use_gaps=False).height
0
off = um.ungap(dr_gapped)
print(off.isgapped, off.gaploc, off.mingap, off.maxgap)
print(um.add_gap(off, off.gaploc, off.mingap, off.maxgap) == dr_gapped)
gone = um.ungap(dr_gapped, delete=True)
print(gone.isgapped, gone.gaploc)
False [6] [1] [5] True False []
Correcting P-values for gaps¶
There is a catch to scanning gapped motifs, however. Each combination of gap sizes is effectively a different motif, and each of them gets its own chance to produce a hit at every position. A motif with five possible gap sizes is thus tested five times as often as one without a gap, and by chance alone can give up to five times as many hits at the same P-value cutoff. The P-values themselves are not wrong, as each one is the P-value of a single match, but they do not account for the number of matches tried.
The correct_gaps argument applies a Bonferroni correction for this: each
P-value is multiplied by the number of gap size combinations and capped at
1, and hits must reach the cutoff after correction. For our motif with five
gap sizes, the P-values of the planted sites are five times larger, though
still far below the cutoff.
um.scan_sequences(dr_gapped, seqs, pvalue=1e-5, correct_gaps=True).select(cols)
| sequence | start | end | strand | score | match | pvalue |
|---|---|---|---|---|---|---|
| str | i64 | i64 | str | f64 | str | f64 |
| "seq1" | 90 | 103 | "+" | 23.364 | "AGGTCA.AGGTCA" | 2.9802e-7 |
| "seq2" | 90 | 104 | "+" | 23.364 | "AGGTCA..AGGTCA" | 2.9802e-7 |
| "seq3" | 90 | 106 | "+" | 23.364 | "AGGTCA....AGGTCA" | 2.9802e-7 |
| "seq4" | 90 | 105 | "-" | 23.364 | "AGGTCA...AGGTCA" | 2.9802e-7 |
For strong sites like these the correction changes little. Where it does matter is for chance hits in long sequences, particularly when gaps are allowed to vary a lot. To see this, let us scan the 100 background sequences from the example data, which contain no planted direct repeats, with the ungapped motif and with two gapped versions of it: one with the gap varying from one to five positions, and another from zero to twenty. Every hit in these sequences is a false positive:
controls = um.example_sequences("controls")
rows = []
for label, motif in [
("no gap", dr),
("gap of 1-5", um.add_gap(dr, gaploc=6, mingap=1, maxgap=5)),
("gap of 0-20", um.add_gap(dr, gaploc=6, mingap=0, maxgap=20)),
]:
rows.append({
"motif": label,
"hits": um.scan_sequences(motif, controls, pvalue=1e-4).height,
"hits_corrected": um.scan_sequences(motif, controls, pvalue=1e-4, correct_gaps=True).height,
})
pl.DataFrame(rows)
| motif | hits | hits_corrected |
|---|---|---|
| str | i64 | i64 |
| "no gap" | 4 | 4 |
| "gap of 1-5" | 12 | 0 |
| "gap of 0-20" | 80 | 3 |
Without the correction, the motif with a gap of one to five positions finds
three times as many chance hits as the ungapped motif, and the motif whose gap
spans 21 possible sizes finds twenty times as many. With the correction,
almost all of them are gone. (Note that the motif without a gap is
unaffected by correct_gaps.) So whether to use the correction depends on
what the hits are for. If you are looking for the best candidate sites, or
comparing hit counts between gapped motifs with different gap ranges or with
ungapped motifs, then the corrected P-values are the fairer ones. If on the
other hand you only want to know which of a handful of sequences contain a
gapped site, then the uncorrected P-values are fine, as long as you keep in
mind that a wide gap range makes chance hits more likely. Either way it is
best to keep gap ranges no wider than needed.
RNA motifs¶
RNA motifs work much the same as DNA motifs, with U in place of T.
switch_alph() switches motifs between the two
alphabets. For example, here is the G-box motif from the example data as an
RNA motif:
gbox = um.example_motifs()[3]
gbox_rna = um.switch_alph(gbox)
print(gbox_rna.alphabet, gbox_rna.letters, gbox_rna.consensus)
RNA ['A', 'C', 'G', 'U'] NCACGUGN
The alphabet of a motif must match that of the sequences it is scanned
against. Scanning the DNA motif against RNA sequences gives no hits at all,
since U is not a DNA letter and so is never part of a hit, and a warning
points this out. The RNA
motif, on the other hand, scans them as expected. Since RNA is usually
single-stranded, it often makes sense to scan only the forward strand with
rc=False, which here leaves a single hit:
rna = um.create_sequences("RNA", seqnum=20, seqlen=200, seed=1)
with warnings.catch_warnings(record=True) as caught:
print(um.scan_sequences(gbox, rna, pvalue=1e-3).height)
print(caught[0].message)
um.scan_sequences(gbox_rna, rna, pvalue=1e-3, rc=False).select(cols)
0 DNA motifs are scanned against sequences that contain U but not T, so they can only match where T is not needed; convert the motifs with switch_alph() or the sequences first
| sequence | start | end | strand | score | match | pvalue |
|---|---|---|---|---|---|---|
| str | i64 | i64 | str | f64 | str | f64 |
| "seq19" | 32 | 40 | "+" | 10.479 | "UCACGUGU" | 0.000122 |
Amino acid motifs¶
Amino acid motifs use the 20 standard amino acids as their alphabet, which is
given as "AA". Let us make one from a few short aligned protein sites. It is
best to always give the alphabet explicitly when creating amino acid motifs,
since letters such as R, K and N are also DNA ambiguity codes, and the
alphabet may otherwise be guessed wrong:
protein_sites = ["LDEIKR", "LNQLKR", "LDELRK", "LNQIKR", "LDEIKK", "LNELRR"]
toy = um.create_motif(protein_sites, alphabet="AA", name="toy", pseudocount=1)
print(toy)
Motif name: toy
Alphabet: AA
Type: PPM
Total IC: 15.18
Pseudocount: 1
Consensus: LBEJKR
Target sites: 6
L B E J K R
A 0.01 0.01 0.01 0.01 0.01 0.01
C 0.01 0.01 0.01 0.01 0.01 0.01
D 0.01 0.44 0.01 0.01 0.01 0.01
E 0.01 0.01 0.58 0.01 0.01 0.01
F 0.01 0.01 0.01 0.01 0.01 0.01
G 0.01 0.01 0.01 0.01 0.01 0.01
H 0.01 0.01 0.01 0.01 0.01 0.01
I 0.01 0.01 0.01 0.44 0.01 0.01
K 0.01 0.01 0.01 0.01 0.58 0.29
L 0.86 0.01 0.01 0.44 0.01 0.01
M 0.01 0.01 0.01 0.01 0.01 0.01
N 0.01 0.44 0.01 0.01 0.01 0.01
P 0.01 0.01 0.01 0.01 0.01 0.01
Q 0.01 0.01 0.29 0.01 0.01 0.01
R 0.01 0.01 0.01 0.01 0.29 0.58
S 0.01 0.01 0.01 0.01 0.01 0.01
T 0.01 0.01 0.01 0.01 0.01 0.01
V 0.01 0.01 0.01 0.01 0.01 0.01
W 0.01 0.01 0.01 0.01 0.01 0.01
Y 0.01 0.01 0.01 0.01 0.01 0.01
The consensus uses a few ambiguity letters of its own. B stands for D or N
and J for I or L, and so the second and fourth positions, which are evenly
split between those pairs, get those letters. Similarly, Z stands for E or Q.
Otherwise the consensus letter is simply the most common amino acid, or X
when the two most common ones are tied or when no amino acid is common at
all. The consensus letter of a single column can be had with
get_consensus_aa(), with the values given
in alphabetical order of the amino acids.
column = [0.0] * 20
column[toy.letters.index("E")] = 0.5
column[toy.letters.index("Q")] = 0.5
um.get_consensus_aa(column)
'Z'
Scanning protein sequences works just like scanning DNA. Let us make five random protein sequences, plant one of the sites in the second and a lower case variant of it in the fourth, and scan them. Amino acids, like DNA and RNA, are matched regardless of case, and matches are always shown in upper case:
proteins = um.create_sequences("AA", seqnum=5, seqlen=300, seed=2)
proteins[1] = plant(proteins[1], "LDEIKR", at=100)
proteins[3] = plant(proteins[3], "ldqlrr", at=200)
um.scan_sequences(toy, proteins, pvalue=1e-5).select(cols)
| sequence | start | end | strand | score | match | pvalue |
|---|---|---|---|---|---|---|
| str | i64 | i64 | str | f64 | str | f64 |
| "seq2" | 100 | 106 | "+" | 20.953 | "LDEIKR" | 6.2500e-8 |
| "seq4" | 200 | 206 | "+" | 18.989 | "LDQLRR" | 4.3750e-7 |
Amino acid motifs can also be compared with one another using
compare_motifs(), as long as all of the
motifs share the same alphabet. Here a second motif made from similar sites
correlates closely with the first, and a random amino acid motif does not:
toy2 = um.create_motif(["LDEIKR", "LDELKR", "LNEIKR"], alphabet="AA", name="toy2", pseudocount=1)
random_aa = um.create_motif(6, alphabet="AA", seed=1, name="random")
um.compare_motifs([toy, toy2, random_aa]).round(3)
array([[1. , 0.926, 0.161],
[0.926, 1. , 0.222],
[0.161, 0.222, 1. ]])
Finally, view_motifs() draws amino acid logos
with the letters coloured by the chemistry of their side chains:
fig = um.view_motifs([toy, toy2])
Custom alphabets¶
Motifs are not limited to biological sequence letters either. Any set of
letters can serve as an alphabet, simply by giving those letters as the
alphabet. For example, protein secondary structure is often written with
three letters, H for helix, E for strand and C for coil, so let us make a
motif of a helix turning into a strand:
ss_sites = ["HHHHCCEEE", "HHHHCCEEE", "HHHCCCEEE", "HHHHCCCEE", "HHHHCEEEE"]
helix_strand = um.create_motif(ss_sites, alphabet="CEH", name="helix-strand", pseudocount=1)
print(helix_strand)
Motif name: helix-strand
Alphabet: CEH
Type: PPM
Total IC: 7.42
Pseudocount: 1
Target sites: 5
1 2 3 4 5 6 7 8 9
C 0.06 0.06 0.06 0.22 0.89 0.72 0.22 0.06 0.06
E 0.06 0.06 0.06 0.06 0.06 0.22 0.72 0.89 0.89
H 0.89 0.89 0.89 0.72 0.06 0.06 0.06 0.06 0.06
There is no consensus for a custom alphabet, since there are no ambiguity letters to build one from, and so the positions are simply numbered. Scanning works as before, but the letters of custom alphabets are matched exactly, including their case. In the following example the third sequence is written in lower case, and so it gets no hits.
structures = {
"prot1": "CCCHHHHHHHCCEEEECCCC",
"prot2": "CCCHHHHCCEEECCCCCCHHHHCCEEEC",
"prot3": "ccchhhhcceeeccc",
}
um.scan_sequences(helix_strand, structures, pvalue=1e-3).select(cols)
| sequence | start | end | strand | score | match | pvalue |
|---|---|---|---|---|---|---|
| str | i64 | i64 | str | f64 | str | f64 |
| "prot1" | 6 | 15 | "+" | 11.835 | "HHHHCCEEE" | 0.000051 |
| "prot1" | 7 | 16 | "+" | 8.437 | "HHHCCEEEE" | 0.000356 |
| "prot2" | 3 | 12 | "+" | 11.835 | "HHHHCCEEE" | 0.000051 |
| "prot2" | 18 | 27 | "+" | 11.835 | "HHHHCCEEE" | 0.000051 |
Logos of custom alphabets are drawn with a categorical palette by default,
which can be changed with the colour_scheme argument:
fig = um.view_motifs(helix_strand, colour_scheme={"H": "#d95f02", "E": "#1b9e77", "C": "#7570b3"})
Strands¶
Only DNA and RNA have a reverse strand. When scanning amino acid or custom
alphabet motifs the rc argument is ignored, which is why all of the hits
above are on the + strand, and similarly compare_motifs()
and view_motifs() only try reverse
complements of DNA and RNA motifs. Asking for the reverse complement of any
other motif with motif_rc() is an error:
try:
um.motif_rc(toy)
except ValueError as e:
print(e)
only DNA/RNA motifs can be reverse complemented
Where next¶
The Scanning large genomes guide shows how to scan genome-sized inputs efficiently, and the API reference describes every function used here in more detail.