Working with motifs¶
This guide introduces the Motif class and the
functions for creating, reading, writing and modifying motifs, and for drawing
them as logos. For an introduction to sequence motifs themselves and the
different ways of representing them, see the
Sequence motifs guide first.
Most of the examples use the twelve synthetic motifs that come with the package, which are built from well-known sites such as the TATA box, the AP-1 site and the plant W-box and G-box.
import universalmotif as um
motifs = um.example_motifs()
[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']
ap1 = motifs[1]
print(ap1)
Motif name: AP-1
Family: bZIP
Organism: synthetic
Alphabet: DNA
Type: PCM
Strands: +-
Total IC: 9.65
Pseudocount: 1
Consensus: NTGASTCAN
Target sites: 49
N T G A S T C A N
A 10 3 0 40 0 2 1 48 16
C 11 1 2 5 24 2 48 1 4
G 9 0 45 2 24 2 0 0 20
T 19 45 2 2 1 43 0 0 9
Only the metadata that is set is shown. All of it is available as attributes:
name: the motif name.altname: an optional alternative name.family: an optional transcription factor or matrix family.organism: an optional organism of origin.matrix: the motif matrix, as a read-only numpy array of shape (letters, positions).alphabetandletters: the alphabet, which is"DNA","RNA","AA"or the letters of a custom alphabet, and its letters in row order.type: one of PCM, PPM, PWM, ICM or CWM.icscore: the total information content of the motif, calculated automatically.nsites: the number of sites the motif was made from, if known.pseudocount: the pseudocount used when converting between types.bkg: the background letter frequencies, as a dict.bkgsites: an optional number of background sites.consensus: the consensus string of DNA, RNA and amino acid motifs, calculated automatically.strand: the strands the motif can be found on, one of"+","-"or"+-".pval,qvalandeval: optional P-, Q- and E-values, such as those from motif discovery.extrainfo: a dict of any other information as strings.width: the number of positions, which is also whatlen()gives.
print(ap1.name, ap1.family, ap1.type, ap1.width, ap1.nsites)
print(ap1.consensus, round(ap1.icscore, 2))
print(ap1.bkg)
ap1.matrix
AP-1 bZIP PCM 9 49.0
NTGASTCAN 9.65
{'A': 0.25, 'C': 0.25, 'G': 0.25, 'T': 0.25}
array([[10., 3., 0., 40., 0., 2., 1., 48., 16.],
[11., 1., 2., 5., 24., 2., 48., 1., 4.],
[ 9., 0., 45., 2., 24., 2., 0., 0., 20.],
[19., 45., 2., 2., 1., 43., 0., 0., 9.]])
Changing metadata¶
Motifs are immutable, so their attributes cannot be assigned to and their matrices cannot be changed in place.
try:
ap1.name = "Jun/Fos"
except AttributeError as e:
print(e)
attribute 'name' of 'universalmotif.Motif' objects is not writable
Instead, replace() returns a copy of the motif with some of its metadata
changed, leaving the original as it was. Any of the metadata can be replaced
this way, but not the matrix, the type, or the consensus and icscore which
are calculated from it:
ap1_new = ap1.replace(altname="Jun/Fos", pval=1e-5,
extrainfo={"source": "example data"})
print(ap1_new)
Motif name: AP-1
Alternate name: Jun/Fos
Family: bZIP
Organism: synthetic
Alphabet: DNA
Type: PCM
Strands: +-
Total IC: 9.65
Pseudocount: 1
Consensus: NTGASTCAN
Target sites: 49
P-value: 1e-5
Extra info: [source] example data
N T G A S T C A N
A 10 3 0 40 0 2 1 48 16
C 11 1 2 5 24 2 48 1 4
G 9 0 45 2 24 2 0 0 20
T 19 45 2 2 1 43 0 0 9
The new values are checked before the copy is made, so that motifs cannot end up with metadata that makes no sense. For example, the strand must be one of the three allowed values and the background must add up to 1.
for change in ({"strand": "x"}, {"bkg": [0.3, 0.2, 0.2, 0.2]}):
try:
ap1.replace(**change)
except ValueError as e:
print(e)
incorrect strand 'x' (must be one of +, -, +-) invalid motif: bkg probabilities must add up to 1
Creating motifs¶
create_motif() makes motifs from consensus
strings, from aligned sites or from matrices, and can also make random ones.
Whatever the input, the type of the returned motif is set with type, which
is a PPM by default, and any of the metadata can be given as well.
From a consensus string¶
A single string is turned into a motif using the IUPAC ambiguity codes, which is handy for making a quick motif to use elsewhere. Let us make one with a pseudocount, and say that it represents 40 sites:
print(um.create_motif("CCNSNGG", nsites=40, pseudocount=1, name="CCNSNGG"))
Motif name: CCNSNGG
Alphabet: DNA
Type: PPM
Strands: +-
Total IC: 8.26
Pseudocount: 1
Consensus: CCNSNGG
Target sites: 40
C C N S N G G
A 0.01 0.01 0.25 0.01 0.25 0.01 0.01
C 0.98 0.98 0.25 0.49 0.25 0.01 0.01
G 0.01 0.01 0.25 0.49 0.25 0.98 0.98
T 0.01 0.01 0.25 0.01 0.25 0.01 0.01
Note that the pseudocount has been applied to the probabilities, which is why the C at position 1 has a probability of 0.98 rather than 1. Without one, the letters that do not match the consensus would have a probability of zero.
From sites¶
A list of strings of the same length is taken as a set of aligned sites,
which are counted, and nsites is set to the number of sites. Let us draw
twenty sites from the W-box motif with
sample_sites() and build a new motif from
them:
wbox = motifs[5]
sites = um.sample_sites(wbox, n=20, seed=1)
print(sites[:5])
print(um.create_motif(sites, type="PCM", name="W-box sites"))
['TTTTGACCG', 'CGTTGACCC', 'ACTTGACTG', 'GACTGACCC', 'CTTTGTCCC']
Motif name: W-box sites
Alphabet: DNA
Type: PCM
Strands: +-
Total IC: 10.04
Pseudocount: 0
Consensus: NNTTGACCC
Target sites: 20
N N T T G A C C C
A 4 2 0 1 0 17 0 0 1
C 7 4 1 0 0 0 18 14 14
G 6 6 0 1 20 0 1 1 5
T 3 8 19 18 0 3 1 5 0
From a matrix¶
A matrix can be given as a numpy array, or anything that can be turned into
one, with the letters as rows and the positions as columns. Its type is
detected from the values. The rows are taken in alphabet order when
alphabet is given, which for RNA is A, C, G and U:
import numpy as np
probabilities = np.array([
[0.7, 0.7, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
[0.1, 0.1, 0.7, 0.7, 0.1, 0.1, 0.1, 0.1],
[0.1, 0.1, 0.1, 0.1, 0.7, 0.7, 0.1, 0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.7, 0.7],
])
print(um.create_motif(probabilities, alphabet="RNA", name="My motif",
pseudocount=1, nsites=20, strand="+"))
Motif name: My motif
Alphabet: RNA
Type: PPM
Strands: +
Total IC: 4.68
Pseudocount: 1
Consensus: AACCGGUU
Target sites: 20
A A C C G G U U
A 0.7 0.7 0.1 0.1 0.1 0.1 0.1 0.1
C 0.1 0.1 0.7 0.7 0.1 0.1 0.1 0.1
G 0.1 0.1 0.1 0.1 0.7 0.7 0.1 0.1
U 0.1 0.1 0.1 0.1 0.1 0.1 0.7 0.7
The type, total IC and consensus are filled in for us. If the rows come in
some other order, give their letters with letters instead and they will be
put into alphabet order.
counts = [[0, 9, 1], # T
[1, 0, 8], # G
[9, 0, 0], # C
[0, 1, 1]] # A
um.create_motif(counts, letters=["T", "G", "C", "A"], type="PCM").matrix
array([[0., 1., 1.],
[9., 0., 0.],
[1., 0., 8.],
[0., 9., 1.]])
Random motifs¶
Without an input, or with just a width, a random motif is made. Each column
is drawn from a Dirichlet distribution parameterized by the background, so it
is worth fixing seed whenever the motif needs to be the same each time:
print(um.create_motif(8, seed=1))
Motif name: motif
Alphabet: DNA
Type: PPM
Strands: +-
Total IC: 9.93
Pseudocount: 0
Consensus: GATCCGWT
Target sites: 189
G A T C C G W T
A 0.07 0.70 0.14 0.01 0.09 0.10 0.30 0.22
C 0.02 0.14 0.01 0.94 0.90 0.01 0.00 0.00
G 0.91 0.16 0.01 0.00 0.00 0.90 0.22 0.00
T 0.00 0.00 0.84 0.05 0.00 0.00 0.48 0.78
Giving a background changes the letters the random motif tends to favour, and the background is kept as the motif's own. Averaged over the positions of a longer motif made with a background rich in C, about 43% of the probability goes to C:
c_rich = um.create_motif(50, bkg={"A": 0.2, "C": 0.4, "G": 0.2, "T": 0.2}, seed=1)
print(c_rich.bkg)
c_rich.matrix.mean(axis=1).round(2)
{'A': 0.2, 'C': 0.4, 'G': 0.2, 'T': 0.2}
array([0.19, 0.43, 0.19, 0.19])
Motifs are not limited to DNA, RNA and amino acids, as alphabet can also be
any set of letters, which are then sorted. Custom alphabets work just as well
with most functions, though they do not have consensus strings or strands:
custom = um.create_motif(6, alphabet="QWERTY", seed=1)
print(custom.alphabet, custom.letters, custom.consensus)
fig = um.view_motifs(custom)
EQRTWY ['E', 'Q', 'R', 'T', 'W', 'Y'] None
Reading and writing motifs¶
There are readers and writers for the MEME, JASPAR, HOMER, TRANSFAC, CIS-BP and UniPROBE formats, and for motifs stored as plain matrices of numbers:
| Format | Reader | Writer |
|---|---|---|
| MEME | read_meme() |
write_meme() |
| JASPAR | read_jaspar() |
write_jaspar() |
| HOMER | read_homer() |
write_homer() |
| TRANSFAC | read_transfac() |
write_transfac() |
| CIS-BP | read_cisbp() |
write_cisbp() |
| UniPROBE | read_uniprobe() |
write_uniprobe() |
| Matrix | read_matrix() |
write_matrix() |
Readers always return a list of motifs, even when a file holds just one. There
is an example file for each format, holding the first four example motifs,
whose path is given by example_file(). Let us
read every one of them and see what comes back:
import polars as pl
readers = {
"meme": um.read_meme,
"jaspar": um.read_jaspar,
"homer": um.read_homer,
"transfac": um.read_transfac,
"cisbp": um.read_cisbp,
"uniprobe": um.read_uniprobe,
"matrix": um.read_matrix,
}
rows = []
for fmt, reader in readers.items():
read = reader(um.example_file(fmt))
rows.append({"format": fmt, "motifs": len(read), "type": read[0].type,
"nsites": read[0].nsites, "family": read[0].family})
pl.DataFrame(rows)
| format | motifs | type | nsites | family |
|---|---|---|---|---|
| str | i64 | str | f64 | str |
| "meme" | 4 | "PPM" | 23.0 | null |
| "jaspar" | 4 | "PCM" | 23.0 | null |
| "homer" | 4 | "PPM" | null | null |
| "transfac" | 4 | "PCM" | 23.0 | "TBP" |
| "cisbp" | 4 | "PPM" | null | "TBP" |
| "uniprobe" | 4 | "PPM" | null | null |
| "matrix" | 4 | "PCM" | 23.0 | null |
Each format stores motifs a little differently, and the readers keep as much
as each one holds. JASPAR, TRANSFAC and plain matrix files hold counts, which
are read as PCMs, while the others hold probabilities and are read as PPMs.
HOMER, CIS-BP and UniPROBE files do not record the number of sites, and only
TRANSFAC and CIS-BP have a field for the family. Anything else worth keeping
goes into extrainfo, such as the score threshold in each HOMER motif header:
um.read_homer(um.example_file("homer"))[0].extrainfo
{'logodds.threshold': '7.46064783005494'}
Writing is the same in reverse. Let us write all twelve example motifs in every format to a temporary directory.
import pathlib
import tempfile
outdir = pathlib.Path(tempfile.mkdtemp())
writers = {
"motifs.meme": um.write_meme,
"motifs.jaspar": um.write_jaspar,
"motifs.homer": um.write_homer,
"motifs.transfac": um.write_transfac,
"motifs.cisbp": um.write_cisbp,
"motifs.uniprobe": um.write_uniprobe,
"motifs.txt": um.write_matrix,
}
for filename, writer in writers.items():
writer(motifs, outdir / filename)
sorted(p.name for p in outdir.iterdir())
['motifs.cisbp', 'motifs.homer', 'motifs.jaspar', 'motifs.meme', 'motifs.transfac', 'motifs.txt', 'motifs.uniprobe']
For example, the start of the JASPAR file looks like this:
print("".join((outdir / "motifs.jaspar").read_text().splitlines(keepends=True)[:10]))
>TATA-box A [ 2 2 23 1 18 11 20 8 10 7 ] C [ 1 0 0 1 3 1 3 1 0 3 ] G [ 11 0 0 1 1 1 0 1 13 5 ] T [ 9 21 0 20 1 10 0 13 0 8 ] >AP-1 A [ 10 3 0 40 0 2 1 48 16 ] C [ 11 1 2 5 24 2 48 1 4 ] G [ 9 0 45 2 24 2 0 0 20 ]
Writers will not overwrite an existing file unless told to with
overwrite=True, and append=True adds the motifs to the end of the file
instead. When appending to a MEME file only the motifs are written, so that
the file keeps a single header:
try:
um.write_meme(motifs, outdir / "motifs.meme")
except FileExistsError:
print("motifs.meme already exists")
um.write_meme(um.create_motif("TTGACY", name="extra"), outdir / "motifs.meme", append=True)
len(um.read_meme(outdir / "motifs.meme"))
motifs.meme already exists
13
Do keep the limitations of each format in mind. For example, all the motifs in a MEME file share one alphabet and one set of background frequencies, and motifs read from files without site numbers will be treated as having 100 sites whenever counts are needed.
Plain matrix files vary the most, so both
write_matrix() and
read_matrix() can be told how the file is laid
out: whether positions are columns or rows, whether the letters are written
too, and how names are marked. Here the motifs are written as probabilities
with one position per row and the names marked with >:
um.write_matrix(motifs[:2], outdir / "rows.txt", type="PPM", positions="rows",
rownames=True, headers=">", digits=3)
print((outdir / "rows.txt").read_text()[:150])
[m.name for m in um.read_matrix(outdir / "rows.txt", positions="rows",
rownames=True, headers=">")]
>TATA-box A C G T 0.094 0.052 0.469 0.385 0.094 0.010 0.010 0.885 0.969 0.010 0.010 0.010 0.052 0.052 0.052 0.844 0.760 0.135 0.052 0.052 0.469 0.052
['TATA-box', 'AP-1']
Converting between types¶
convert_type() turns motifs into PCMs, PPMs,
PWMs or ICMs. It takes a single motif or a list, as do most functions that
change motifs. The example motifs are counts, so let us convert the CRE motif
to probabilities.
cre = motifs[2]
print(cre.type, cre.pseudocount)
um.convert_type(cre, "PPM").matrix.round(3)
PCM 1.0
array([[0.069, 0.027, 0.92 , 0.069, 0.005, 0.027, 0.005, 0.793],
[0.048, 0.027, 0.005, 0.92 , 0.069, 0.027, 0.835, 0.176],
[0.069, 0.941, 0.048, 0.005, 0.878, 0.027, 0.09 , 0.005],
[0.814, 0.005, 0.027, 0.005, 0.048, 0.92 , 0.069, 0.027]])
The motif's pseudocount of 1 is used for every conversion, which keeps the
PWM finite even where a letter was never seen. A different pseudocount can be
passed to convert_type(), and without one the
zero counts become minus infinity:
print(um.convert_type(cre, "PWM").matrix.round(2))
um.convert_type(cre, "PWM", pseudocount=0).matrix.round(2)
[[-1.85 -3.23 1.88 -1.85 -5.55 -3.23 -5.55 1.66] [-2.38 -3.23 -5.55 1.88 -1.85 -3.23 1.74 -0.51] [-1.85 1.91 -2.38 -5.55 1.81 -3.23 -1.47 -5.55] [ 1.7 -5.55 -3.23 -5.55 -2.38 1.88 -1.85 -3.23]]
array([[-1.94, -3.52, 1.9 , -1.94, -inf, -3.52, -inf, 1.69],
[-2.52, -3.52, -inf, 1.9 , -1.94, -3.52, 1.76, -0.52],
[-1.94, 1.94, -2.52, -inf, 1.83, -3.52, -1.52, -inf],
[ 1.72, -inf, -3.52, -inf, -2.52, 1.9 , -1.94, -3.52]])
Converting to counts needs the number of sites the motif was made from. When
nsites is not known, as with the motifs read from HOMER files, 100 sites are
assumed, so do keep in mind that these counts are not a measured sample size:
homer_motif = um.read_homer(um.example_file("homer"))[0]
print(homer_motif.nsites)
um.convert_type(homer_motif, "PCM").matrix.sum(axis=0)
None
array([100., 100., 100., 100., 100., 100., 100., 100., 100., 100.])
Note that changing a motif's pseudocount with replace() does not itself
change its matrix, since the pseudocount is only applied when the motif is
converted to another type. normalize() applies
it straight away. Converting to an ICM can also use relative entropy, which is
explained in the Sequence motifs guide along with the
other types.
Summarizing and filtering motifs¶
summarize_motifs() gives an overview of a
list of motifs as a polars DataFrame, with one row per motif and the columns
that are empty for every motif left out:
summary = um.summarize_motifs(motifs)
with pl.Config(tbl_rows=12):
display(summary)
| name | family | organism | consensus | alphabet | strand | icscore | nsites |
|---|---|---|---|---|---|---|---|
| str | str | str | str | str | str | f64 | f64 |
| "TATA-box" | "TBP" | "synthetic" | "KTATAWAWRN" | "DNA" | "+-" | 8.82611 | 23.0 |
| "AP-1" | "bZIP" | "synthetic" | "NTGASTCAN" | "DNA" | "+-" | 9.646674 | 49.0 |
| "CRE" | "bZIP" | "synthetic" | "TGACGTCA" | "DNA" | "+-" | 10.72504 | 46.0 |
| "G-box" | "bZIP" | "synthetic" | "NCACGTGN" | "DNA" | "+-" | 7.899263 | 71.0 |
| "E-box" | "bHLH" | "synthetic" | "NCAGCTGN" | "DNA" | "+-" | 8.214013 | 33.0 |
| "W-box" | "WRKY" | "synthetic" | "GBTTGACYS" | "DNA" | "+-" | 9.797764 | 78.0 |
| "GATA" | "GATA" | "synthetic" | "CWGATART" | "DNA" | "+-" | 6.740901 | 53.0 |
| "GC-box" | "C2H2 zinc finger" | "synthetic" | "GGGGCGGGG" | "DNA" | "+-" | 11.75613 | 43.0 |
| "ETS" | "ETS" | "synthetic" | "NCCGGAAGTG" | "DNA" | "+-" | 10.570188 | 63.0 |
| "CArG-box" | "MADS" | "synthetic" | "CCWWTWWTGG" | "DNA" | "+-" | 8.986613 | 24.0 |
| "Telo-box" | "Myb-related" | "synthetic" | "AAACCCTAA" | "DNA" | "+-" | 11.240669 | 39.0 |
| "CCAAT-box" | "NF-Y" | "synthetic" | "CNCCAATNH" | "DNA" | "+-" | 7.712225 | 44.0 |
filter_motifs() keeps the motifs that pass
every filter it is given. Text filters, such as family or name, keep
motifs with any of the values given, while width, icscore and nsites are
minimums:
print([m.name for m in um.filter_motifs(motifs, family="bZIP")])
[m.name for m in um.filter_motifs(motifs, icscore=9, width=9)]
['AP-1', 'CRE', 'G-box']
['AP-1', 'W-box', 'GC-box', 'ETS', 'Telo-box']
For anything more elaborate, the summary DataFrame can be filtered with polars instead and the motifs then picked out by name. For example, let us find the motifs with an ACGT core in their consensus:
acgt = summary.filter(pl.col("consensus").str.contains("ACGT"))["name"]
[m.name for m in um.filter_motifs(motifs, name=acgt)]
['CRE', 'G-box']
Merging motifs¶
merge_motifs() combines several motifs into
one. The motifs are first aligned to the one with the highest information
content, in the same way as compare_motifs()
aligns them, and the aligned probabilities are then averaged. Let us make
three overlapping motifs and draw them.
m1 = um.create_motif("TTAAACCCC", name="1")
m2 = um.create_motif("AACC", name="2")
m3 = um.create_motif("AACCCCGG", name="3")
fig = um.view_motifs([m1, m2, m3], show_positions_once=False)
Motifs that cannot overlap the others by at least min_overlap positions are
left out of the merge with a warning, which with the default of 5 positions is
the case for the second motif. The name of the merged motif tells us which
motifs went into it:
import warnings
with warnings.catch_warnings(record=True) as caught:
print(um.merge_motifs([m1, m2, m3]).name)
print(caught[0].message)
merged = um.merge_motifs([m1, m2, m3], min_overlap=4)
print(merged.name)
fig = um.view_motifs(merged)
merged_1+3 1 motif could not overlap the anchor by 5 columns and was left out: 2 merged_1+2+3
For larger sets of motifs, merge_similar()
finds groups of similar motifs and merges each group, which is useful for
reducing redundancy, for example when the same motif comes from several
databases. Let us mimic this by building a second version of each example
motif from thirty sites sampled from it. With return_clusters=True we get
the groups instead of the merged motifs:
resampled = [
um.create_motif(um.sample_sites(m, n=30, seed=i), type="PCM", pseudocount=1,
name=m.name + " (resampled)")
for i, m in enumerate(motifs)
]
clusters = um.merge_similar(motifs + resampled, return_clusters=True)
clusters.select("name", "cluster").head(4).vstack(
clusters.select("name", "cluster").slice(12, 4))
| name | cluster |
|---|---|
| str | i64 |
| "TATA-box" | 0 |
| "AP-1" | 1 |
| "CRE" | 2 |
| "G-box" | 3 |
| "TATA-box (resampled)" | 0 |
| "AP-1 (resampled)" | 1 |
| "CRE (resampled)" | 2 |
| "G-box (resampled)" | 3 |
Each motif has been grouped with its resampled version and with nothing else, so merging gives back twelve motifs.
merged = um.merge_similar(motifs + resampled)
print(len(merged))
merged[0].name
12
'merged_TATA-box+TATA-box (resampled)'
Do keep in mind that merging motifs which are not actually alike gives a motif that represents neither of them well. The AP-1 and CRE sites are related, for example, but merging them gives a rather muddled consensus:
um.merge_motifs([motifs[1], motifs[2]]).consensus
'TKRMSTCAN'
Reverse complements¶
motif_rc() gives the reverse complement of DNA
and RNA motifs, as for the TATA-box motif here:
tata = motifs[0]
print(tata.consensus)
print(um.motif_rc(tata).consensus)
KTATAWAWRN NYWTWTATAM
Many sites are palindromes, so that their reverse complements are the same as the originals. The consensus of the AP-1 motif reads the same either way, for example, though its matrix still changes slightly since the two halves of the example motif are not identical:
print(um.motif_rc(ap1).consensus)
um.motif_rc(ap1).matrix
NTGASTCAN
array([[ 9., 0., 0., 43., 1., 2., 2., 45., 19.],
[20., 0., 0., 2., 24., 2., 45., 0., 9.],
[ 4., 1., 48., 2., 24., 5., 2., 1., 11.],
[16., 48., 1., 2., 0., 40., 0., 3., 10.]])
Trimming motifs¶
Motifs often have uninformative positions at their edges.
trim_motifs() removes positions from the edges
until it reaches one with an information content of at least min_ic, which
is 0.25 bits by default. Trimming can also be limited to one side with
trim_from:
loose = um.create_motif("NNGCSGCGGNN")
print(um.trim_motifs(loose).consensus)
print(um.trim_motifs(loose, trim_from="right").consensus)
GCSGCGG NNGCSGCGG
Several of the example motifs have weak edges, which raising min_ic trims
further still.
with pl.Config(tbl_rows=12):
display(pl.DataFrame({
"name": [m.name for m in motifs],
"original": [m.consensus for m in motifs],
"trimmed": [m.consensus for m in um.trim_motifs(motifs)],
"min_ic=0.5": [m.consensus for m in um.trim_motifs(motifs, min_ic=0.5)],
}))
| name | original | trimmed | min_ic=0.5 |
|---|---|---|---|
| str | str | str | str |
| "TATA-box" | "KTATAWAWRN" | "KTATAWAWR" | "TATAWAWR" |
| "AP-1" | "NTGASTCAN" | "TGASTCA" | "TGASTCA" |
| "CRE" | "TGACGTCA" | "TGACGTCA" | "TGACGTCA" |
| "G-box" | "NCACGTGN" | "CACGTG" | "CACGTG" |
| "E-box" | "NCAGCTGN" | "CAGCTG" | "CAGCTG" |
| "W-box" | "GBTTGACYS" | "GBTTGACYS" | "TTGACYS" |
| "GATA" | "CWGATART" | "CWGATART" | "GATA" |
| "GC-box" | "GGGGCGGGG" | "GGGGCGGGG" | "GGGGCGGGG" |
| "ETS" | "NCCGGAAGTG" | "CCGGAAGT" | "CCGGAAGT" |
| "CArG-box" | "CCWWTWWTGG" | "CCWWTWWTGG" | "CCWWTWWTGG" |
| "Telo-box" | "AAACCCTAA" | "AAACCCTAA" | "AAACCCTAA" |
| "CCAAT-box" | "CNCCAATNH" | "CNCCAAT" | "CCAAT" |
A motif with no position above min_ic is trimmed away entirely, which is an
error for a single motif; in a list it is dropped with a warning instead.
Contribution weight matrices have their own trimming function,
trim_cwm().
Rounding motifs¶
round_motif() sets small probabilities to zero
and renormalizes each position, which can tidy up motifs with a lot of
near-zero noise. Probabilities below 5% are rounded off by default. Let us
make a motif with a large pseudocount and round it:
noisy = um.create_motif("ATCGATGC", pseudocount=10, nsites=100, name="noisy")
rounded = um.round_motif(noisy).replace(name="rounded")
print(noisy.matrix.round(3)[:, :2])
print(rounded.matrix[:, :2])
fig = um.view_motifs([noisy, rounded], use_type="PPM", rc=False)
[[0.932 0.023] [0.023 0.023] [0.023 0.023] [0.023 0.932]] [[1. 0.] [0. 0.] [0. 0.] [0. 1.]]
Switching between DNA and RNA¶
Since not every format or tool supports RNA motifs,
switch_alph() switches motifs between the DNA
and RNA alphabets, in either direction.
tata_rna = um.switch_alph(tata)
print(tata_rna.alphabet, tata_rna.letters, tata_rna.consensus)
um.switch_alph(tata_rna).consensus
RNA ['A', 'C', 'G', 'U'] KUAUAWAWRN
'KTATAWAWRN'
Drawing logos¶
view_motifs() draws motifs as sequence logos
with matplotlib. By default motifs are drawn as ICMs, with the height of each
position given by its information content in bits:
fig = um.view_motifs(wbox)
Several motifs are aligned to the one with the highest information content before being drawn one above the other, and any that align better as reverse complements are flipped, with " [RC]" added to their names. The names can be put beside the logos instead of above them, and each logo can get its own position numbers.
bzip = um.filter_motifs(motifs, family="bZIP")
fig = um.view_motifs(bzip, names_pos="right", show_positions_once=False)
The G-box has been reverse complemented here even though it is a palindrome,
since the flipped version happens to line up slightly better with the others.
rc=False turns this off, and sort_by changes the order of the logos.
Logos can be drawn as any type with use_type, as the
Sequence motifs guide shows. Since
view_motifs() returns a matplotlib figure, it
can be adjusted afterwards like any other figure before being saved, for
example by rotating the position numbers of a long motif:
fig = um.view_motifs(um.create_motif(15, seed=2, name="random"))
for ax in fig.axes:
ax.tick_params(axis="x", labelrotation=90)
Finally, view_logo() draws a logo from any
numeric matrix, which makes it useful for all sorts of text logos. The heights
can be anything, including negative values, and each row can be named with
more than one letter. Let us take the CRE motif as an ICM, change a couple of
heights, and keep just two of its rows.
heights = um.convert_type(cre, "ICM").matrix.copy()
heights[0, 3] = 2 # A at position 4
heights[3, 4] = -0.5 # T at position 5
heights = heights[[0, 3], :]
heights.round(2)
array([[ 0.07, 0.04, 1.38, 2. , 0.01, 0.04, 0.01, 0.88],
[ 0.83, 0.01, 0.04, 0.01, -0.5 , 1.35, 0.08, 0.03]])
ax = um.view_logo(heights, letters=["AA", "T"],
colour_scheme={"A": "#109648", "T": "#D62839"})
Note that view_logo() draws onto a set of
matplotlib axes and returns them, so that several logos can be placed in one
figure by passing ax.
Where next¶
The API reference describes everything covered here in full, along with the functions for scanning sequences, comparing motifs and testing for enrichment.