Getting started¶
This guide is a quick tour of the package: what a motif looks like here, how to convert between its different forms, how to read and write motif files, and how to use motifs to scan and compare. Each section only scratches the surface, and the API reference has the details.
Everything is imported from the top-level package:
import universalmotif as um
Motifs¶
A motif can be made from a consensus string, from a set of aligned sites, or from a matrix. For example, here is one made from three sites:
sites = ["TGACTCA", "TGAGTCA", "TGACTAA"]
ap1 = um.create_motif(sites, pseudocount=1, name="AP-1")
print(ap1)
Motif name: AP-1
Alphabet: DNA
Type: PPM
Strands: +-
Total IC: 6.05
Pseudocount: 1
Consensus: TGASTMA
Target sites: 3
T G A S T M A
A 0.06 0.06 0.81 0.06 0.06 0.31 0.81
C 0.06 0.06 0.06 0.56 0.06 0.56 0.06
G 0.06 0.81 0.06 0.31 0.06 0.06 0.06
T 0.81 0.06 0.06 0.06 0.81 0.06 0.06
Printing a motif shows its metadata and its matrix, with one column per position and the consensus letters as column names. Note that the consensus uses the IUPAC ambiguity codes, so the S at position 4 means C or G. The metadata and matrix are all available as attributes:
print(ap1.name, ap1.alphabet, ap1.type, ap1.consensus)
print(round(ap1.icscore, 2))
ap1.matrix.round(2)
AP-1 DNA PPM TGASTMA 6.05
array([[0.06, 0.06, 0.81, 0.06, 0.06, 0.31, 0.81],
[0.06, 0.06, 0.06, 0.56, 0.06, 0.56, 0.06],
[0.06, 0.81, 0.06, 0.31, 0.06, 0.06, 0.06],
[0.81, 0.06, 0.06, 0.06, 0.81, 0.06, 0.06]])
In a Jupyter notebook, a motif on its own at the end of a cell is shown
as a logo next to its metadata, with the matrix tucked away underneath.
logo_svg() gives the same logo as an SVG
image for reports and web pages:
ap1
| Name | AP-1 |
|---|---|
| Alphabet | DNA |
| Type | PPM |
| Strands | +- |
| Total IC | 6.05 |
| Pseudocount | 1 |
| Consensus | TGASTMA |
| Target sites | 3 |
Matrix
T G A S T M A A 0.06 0.06 0.81 0.06 0.06 0.31 0.81 C 0.06 0.06 0.06 0.56 0.06 0.56 0.06 G 0.06 0.81 0.06 0.31 0.06 0.06 0.06 T 0.81 0.06 0.06 0.06 0.81 0.06 0.06
Motifs are immutable, so functions always return new motifs rather than
changing the ones they are given. To change some of the metadata, use
replace():
ap1 = ap1.replace(family="bZIP", organism="Homo sapiens")
print(ap1.family, ap1.organism)
bZIP Homo sapiens
Motif types¶
The same motif can be represented as counts, probabilities, log-odds
weights or information content, which are the PCM, PPM, PWM and ICM types
respectively. Converting between them
is a single call to convert_type(). Converting to counts uses the number
of sites the motif was made from, and converting to weights uses its
pseudocount to keep zero probabilities from becoming minus infinity:
pwm = um.convert_type(ap1, "PWM")
pwm.matrix.round(2)
array([[-2. , -2. , 1.7 , -2. , -2. , 0.32, 1.7 ],
[-2. , -2. , -2. , 1.17, -2. , 1.17, -2. ],
[-2. , 1.7 , -2. , 0.32, -2. , -2. , -2. ],
[ 1.7 , -2. , -2. , -2. , 1.7 , -2. , -2. ]])
um.convert_type(ap1, "ICM").matrix.round(2)
array([[0.06, 0.06, 0.82, 0.03, 0.06, 0.16, 0.82],
[0.06, 0.06, 0.06, 0.29, 0.06, 0.29, 0.06],
[0.06, 0.82, 0.06, 0.16, 0.06, 0.03, 0.06],
[0.82, 0.06, 0.06, 0.03, 0.82, 0.03, 0.06]])
Reading and writing motifs¶
There are readers and writers for the MEME, JASPAR, HOMER, TRANSFAC, CIS-BP, UniPROBE and plain matrix formats. Let us write two motifs to a MEME file and read them back. Readers always return a list, even for a single motif:
import tempfile, pathlib
creb_sites = ["TGACGTCA", "TGACGTCA", "TTACGTAA"]
creb = um.create_motif(creb_sites, pseudocount=1, name="CREB")
path = pathlib.Path(tempfile.mkdtemp()) / "motifs.meme"
um.write_meme([ap1, creb], path)
motifs = um.read_meme(path)
[m.name for m in motifs]
['AP-1', 'CREB']
summarize_motifs() gives an overview of a list of motifs as a
DataFrame, which is handy once there are more than a few of them:
um.summarize_motifs(motifs)
| name | consensus | alphabet | strand | icscore | nsites |
|---|---|---|---|---|---|
| str | str | str | str | f64 | f64 |
| "AP-1" | "TGASTMA" | "DNA" | "+-" | 6.050407 | 3.0 |
| "CREB" | "TKACGTMA" | "DNA" | "+-" | 7.057014 | 3.0 |
Changing motifs¶
A few of the functions for working with motifs. motif_rc() reverse
complements a motif, trim_motifs() removes uninformative positions from
its edges, and merge_motifs() aligns several motifs and combines them into
one:
print(um.motif_rc(ap1).consensus)
loose = um.create_motif("NNTGACTCANN", pseudocount=1)
print(um.trim_motifs(loose).consensus)
longer = um.create_motif(["ATGAGTCAT"], pseudocount=1)
print(um.merge_motifs([ap1, longer], min_overlap=5).consensus)
TKASTCA TGACTCA ATGAGTCAT
Logos¶
view_motifs() draws sequence logos with matplotlib, aligning several
motifs to one another first and reverse complementing them where that
lines them up better:
fig = um.view_motifs([ap1, creb])
Scanning sequences¶
Let us make some random sequences and plant an AP-1 site in a few of them,
then scan for both motifs. scan_sequences() returns one row per hit, each
with its exact P-value; only hits at or below the pvalue cutoff are
kept:
import polars as pl
seqs = um.create_sequences(seqnum=5, seqlen=200, seed=1)
for i in (0, 2, 3):
seqs[i] = seqs[i][:100] + "TGACTCA" + seqs[i][107:]
hits = um.scan_sequences([ap1, creb], seqs, pvalue=1e-3)
with pl.Config(tbl_rows=20):
display(hits.select("motif", "sequence", "start", "strand", "match", "pvalue"))
| motif | sequence | start | strand | match | pvalue |
|---|---|---|---|---|---|
| str | str | i64 | str | str | f64 |
| "AP-1" | "seq1" | 100 | "+" | "TGACTCA" | 0.000061 |
| "AP-1" | "seq1" | 100 | "-" | "TGAGTCA" | 0.000183 |
| "AP-1" | "seq1" | 115 | "+" | "TGATTCA" | 0.000488 |
| "AP-1" | "seq1" | 115 | "-" | "TGAATCA" | 0.000488 |
| "AP-1" | "seq3" | 100 | "+" | "TGACTCA" | 0.000061 |
| "AP-1" | "seq3" | 100 | "-" | "TGAGTCA" | 0.000183 |
| "AP-1" | "seq4" | 100 | "+" | "TGACTCA" | 0.000061 |
| "AP-1" | "seq4" | 100 | "-" | "TGAGTCA" | 0.000183 |
| "CREB" | "seq3" | 129 | "+" | "TGAGGTCA" | 0.000397 |
| "CREB" | "seq3" | 129 | "-" | "TGACCTCA" | 0.000397 |
| "CREB" | "seq3" | 192 | "+" | "AGACGTCA" | 0.000397 |
| "CREB" | "seq3" | 192 | "-" | "TGACGTCT" | 0.000397 |
Each planted site is found on both strands, since TGACTCA reads as TGAGTCA on the reverse strand, which matches the motif almost as well. The other hits are chance matches in the random sequence: one AP-1-like site in seq1 and two CREB-like ones in seq3, each of them close enough to palindromic to turn up on both strands too. A handful of these is to be expected at a cutoff as lenient as 1e-3, which is why stricter cutoffs such as 1e-4 or 1e-5 are more usual for real sequences.
Sequences can also be read from a FASTA file with read_fasta(), and very
large ones, such as a whole genome, can be scanned straight from the file
by passing its path to scan_sequences().
Comparing motifs¶
Finally, compare_motifs() measures how similar motifs are to one another
by the Pearson correlation of their aligned columns, giving a matrix by
default or a table of the best alignments with P-values:
like_sites = ["TGAGTCA", "TGAGTCA", "TGACTCT"]
ap1_like = um.create_motif(like_sites, pseudocount=1, name="AP-1 like")
um.compare_motifs([ap1, ap1_like, creb]).round(3)
array([[1. , 0.963, 0.437],
[0.963, 1. , 0.446],
[0.437, 0.446, 1. ]])
Where next¶
This only covers the basics, and the other guides go further. Sequence motifs explains the motif types and how they relate to one another, Working with motifs and Working with sequences cover the everyday tools, and Scanning and enrichment, Comparing motifs and P-values and Choosing a background go into the analyses. For whole genomes, see Scanning large genomes, and A ChIP-seq workflow puts most of this together on real peaks, and Building a motif library gathers motifs from public databases. The API reference describes every function.