Working with sequences¶
This guide goes through the functions for handling sequences: reading and writing FASTA files, generating random sequences, calculating backgrounds, shuffling, measuring sequence complexity, working in sliding windows, masking, and sampling sites from motifs to plant in sequences. Most of these are a means to an end, since scanning and enrichment analyses work best with sequences that have been prepared with some care and with backgrounds that match them. Those analyses are covered in Scanning and enrichment.
There is no special class for sequences. A single sequence is a string, a set of sequences is a list of strings or a dict from names to sequences, and functions which return sequences give them back in the same form they were given.
import universalmotif as um
import polars as pl
Reading and writing FASTA files¶
read_fasta() reads a FASTA file into a dict
from names to sequences, where each name is the whole header line. Gzipped
files are recognized from their first bytes, so there is nothing special to
do for them. The package comes with a gzipped FASTA file of 100 example
sequences, each 1000 bases long:
path = um.example_file("fasta")
seqs = um.read_fasta(path)
print(path.name, len(seqs))
first = next(iter(seqs))
print(first, len(seqs[first]))
print(seqs[first][:60])
example_targets.fa.gz 100 target_001 1000 TAGGGTACAATATATTTCTATCAAAAGACAGGAATCACCTTAACCACCTAATGGGCTGAG
These are the same sequences that
example_sequences() returns, and they
will be used throughout this guide.
Writing sequences is done with
write_fasta(), which wraps lines at 80
characters by default. Let us write the first three sequences with shorter
lines and have a look at the start of the file:
import tempfile, pathlib
tmp = pathlib.Path(tempfile.mkdtemp())
three = dict(list(seqs.items())[:3])
um.write_fasta(three, tmp / "three.fa", width=60)
print("\n".join((tmp / "three.fa").read_text().splitlines()[:4]))
>target_001 TAGGGTACAATATATTTCTATCAAAAGACAGGAATCACCTTAACCACCTAATGGGCTGAG TGTTCTAGGGTAAGTTTTCTAAAGTAATGTTTCGTAAAGTGTTAGTATAGACCTGGCCTA AATTTTTCATTAAATATAAAAATTGGTATAAAGACTAAAGAACTATTGATCCTCTATAGT
Files are always written as plain text, whatever their extension. To get a
gzipped file instead, compress it afterwards with Python's gzip module,
after which it can be read back like any other:
import gzip, shutil
with open(tmp / "three.fa", "rb") as src, gzip.open(tmp / "three.fa.gz", "wb") as dst:
shutil.copyfileobj(src, dst)
um.read_fasta(tmp / "three.fa.gz") == three
True
Creating random sequences¶
create_sequences() generates a number
of random sequences of a given length. By default these are DNA with every
letter equally likely, but the alphabet can also be RNA, amino acids
("AA") or any set of letters you like. Since the results are random,
setting a seed makes them reproducible:
um.create_sequences(seqnum=3, seqlen=50, seed=1)
['TTTGGACAGCGCATCTCGAATGATGTGGAATCGGGGCTTTGTCCAAAGTA', 'TGGTAGCAAGTATTGAGCTGCTTCTTACTCGAGTAAAGTATGGCGACGTG', 'ATAGAAGGACGTCTGTCTAGAGTAGGTAGGGGCAGGCACCTGCACCTTAG']
print(um.create_sequences("AA", seqnum=2, seqlen=50, seed=1))
print(um.create_sequences("QWERTY", seqnum=2, seqlen=50, seed=1))
['WVWPNEHFQHMIETLSGPDFYPCSMYQNCATLRMMMHYWYMYKIADDMVD', 'SNMYDPKDEPSFTSPAMLVRGVTIYVFGWLQCMWFEANVAYMMKNCIPSN'] ['YYYTTEQQWQTREWRWQTEQYTEWTYWTEEWRWTTTQYYYTYRREEETYE', 'WTTYETREQTWQWWTETRYWQYWRYYQQYRWETYQQETYEYTTRTERTWT']
The frequencies of the letters are set with freqs, for example to make
AT-rich sequences:
um.create_sequences(seqnum=3, seqlen=50, seed=1,
freqs={"A": 0.35, "C": 0.15, "G": 0.15, "T": 0.35})
['TTTGGAAATAGCATCTAGAATGATGTTGAATCTGGGATTTGTCCAAAGTA', 'TGGTAGCAAGTATTGAGCTTATTCTTAATCTAGTAAAGTATGGCGACGTG', 'ATATAAGTAATTATTTCTAGAGTATTTATTGGCATGCACCTTCAACTTAG']
freqs can also give the frequency of every dinucleotide, or of every
k-mer for a larger k. In this case the sequences are built in a
Markov-style manner, where the probability of each letter depends on the
letters which precede it. Let us use this to make 100 sequences in which
the dinucleotide CG is rare, much as it is in many genomes, as these will
come in handy once we get to shuffling:
kmers = um.get_kmers("ACGT", 2)
cg_freqs = {kmer: 0.01 if kmer == "CG" else 0.99 / 15 for kmer in kmers}
cg_poor = um.create_sequences(seqnum=100, seqlen=1000, freqs=cg_freqs, seed=1)
targets = um.example_sequences("targets")
um.get_bkg(targets, k=1)
| kmer | count | probability |
|---|---|---|
| str | i64 | f64 |
| "A" | 32793 | 0.32793 |
| "C" | 17095 | 0.17095 |
| "G" | 17101 | 0.17101 |
| "T" | 33011 | 0.33011 |
These sequences are AT-rich, with A and T each making up about a third of the letters. Several k-mer sizes can be counted at once, in which case the rows are sorted by size and then alphabetically. Let us count the mono- and dinucleotides of the CG-poor sequences:
with pl.Config(tbl_rows=20):
display(um.get_bkg(cg_poor, k=[1, 2]))
| kmer | count | probability |
|---|---|---|
| str | i64 | f64 |
| "A" | 26873 | 0.26873 |
| "C" | 26595 | 0.26595 |
| "G" | 19626 | 0.19626 |
| "T" | 26906 | 0.26906 |
| "AA" | 6755 | 0.067618 |
| "AC" | 6694 | 0.067007 |
| "AG" | 6642 | 0.066486 |
| "AT" | 6749 | 0.067558 |
| "CA" | 8444 | 0.084525 |
| "CC" | 8413 | 0.084214 |
| "CG" | 1321 | 0.013223 |
| "CT" | 8395 | 0.084034 |
| "GA" | 4886 | 0.048909 |
| "GC" | 4843 | 0.048478 |
| "GG" | 4902 | 0.049069 |
| "GT" | 4972 | 0.04977 |
| "TA" | 6763 | 0.067698 |
| "TC" | 6624 | 0.066306 |
| "TG" | 6735 | 0.067417 |
| "TT" | 6762 | 0.067688 |
Note how rare CG is compared with the other dinucleotides starting with a C: there are 1321 of them, against more than 8000 each of CA, CC and CT. There are fewer Gs overall as well, since a G is unlikely to follow a C.
get_kmers() lists every k-mer of some letters
in the same order, and count_kmers() counts
the k-mers in a single string of any alphabet. When no alphabet is given,
it is taken from the letters in the string, so k-mers which never occur
still get a row:
print(um.get_kmers("ACGT", 2))
um.count_kmers("DASDSDDSASDSSA", k=2)
['AA', 'AC', 'AG', 'AT', 'CA', 'CC', 'CG', 'CT', 'GA', 'GC', 'GG', 'GT', 'TA', 'TC', 'TG', 'TT']
| kmer | count |
|---|---|
| str | i64 |
| "AA" | 0 |
| "AD" | 0 |
| "AS" | 2 |
| "DA" | 1 |
| "DD" | 1 |
| "DS" | 3 |
| "SA" | 2 |
| "SD" | 3 |
| "SS" | 1 |
A few more options of get_bkg() are worth
knowing about. merge=False gives every sequence its own counts, rc=True
counts the reverse complements as well, and pseudocount is added to every
count. Finally, to_meme writes the background to a file in the format
used by MEME background files instead of returning it, which needs a
pseudocount whenever some k-mer is never seen:
um.get_bkg(targets, k=[1, 2], pseudocount=1, to_meme=tmp / "targets.bkg")
print("\n".join((tmp / "targets.bkg").read_text().splitlines()[:7]))
# order 0 A 3.279e-01 C 1.710e-01 G 1.710e-01 T 3.301e-01 # order 1 AA 1.088e-01
Per-sequence counts are also a simple way to compare sequences by their composition. As an example, let us generate three groups of 20 sequences with different letter frequencies, calculate the trinucleotide frequencies of every sequence, and plot the first two principal components:
import numpy as np
import matplotlib.pyplot as plt
groups = {
"AT-rich": {"A": 0.35, "C": 0.15, "G": 0.15, "T": 0.35},
"AC-rich": {"A": 0.35, "C": 0.35, "G": 0.15, "T": 0.15},
"uniform": None,
}
rows, labels = [], []
for i, (label, freqs) in enumerate(groups.items()):
group = um.create_sequences(seqnum=20, seqlen=500, freqs=freqs, seed=i)
bkg = um.get_bkg(group, k=3, merge=False)
wide = bkg.pivot(on="kmer", index="sequence", values="probability")
rows.append(wide.drop("sequence").to_numpy())
labels += [label] * len(group)
x = np.vstack(rows)
u, d, _ = np.linalg.svd(x - x.mean(axis=0), full_matrices=False)
pcs = u[:, :2] * d[:2]
labels = np.array(labels)
fig, ax = plt.subplots(figsize=(4.5, 4))
for label in groups:
ax.scatter(*pcs[labels == label].T, s=15, label=label)
ax.set_xlabel("PC1")
ax.set_ylabel("PC2")
ax.legend();
The three groups separate cleanly. Real sequences will rarely be this tidy of course, and trinucleotides are only one choice amongst many; the k-mer size is up to you.
Shuffling sequences¶
Motif enrichment analyses and de novo motif discovery compare a set of target sequences against a set of background sequences. For such a comparison to say anything about motifs, the two sets should share the same overall composition, otherwise differences in letter frequencies alone can make motifs look enriched or depleted. When no suitable background sequences are at hand, the usual approach is to shuffle the targets in a way that keeps their k-mer composition, such as their dinucleotide frequencies.
shuffle_sequences() offers three
methods for this. The default, "euler", keeps the counts of every k-mer
exactly, and with them the counts of every letter. The first and last k-1
letters stay where they are though, as a consequence of how the method
works. The "markov" method instead calculates the k-mer frequencies and
generates new sequences from them, much as
create_sequences() does, so the
frequencies are only kept approximately and the letter counts will likely
change a little. Lastly, "linear" splits each sequence into blocks of k
letters and puts them back together in random order. This keeps the letter
counts exactly, but only the k-mers within the blocks are kept, not those
which span two of them.
Let us shuffle the CG-poor sequences in each of these ways and compare the counts of the letters and of a few dinucleotides:
shuffled = {
"euler, k=1": um.shuffle_sequences(cg_poor, k=1, seed=1),
"euler, k=2": um.shuffle_sequences(cg_poor, k=2, seed=1),
"markov, k=2": um.shuffle_sequences(cg_poor, k=2, method="markov", seed=1),
"linear, k=2": um.shuffle_sequences(cg_poor, k=2, method="linear", seed=1),
}
counts = um.get_bkg(cg_poor, k=[1, 2]).select("kmer", original="count")
for label, s in shuffled.items():
counts = counts.with_columns(um.get_bkg(s, k=[1, 2])["count"].alias(label))
counts.filter(pl.col("kmer").is_in(["A", "C", "G", "T", "CA", "CG", "GC", "TG"]))
| kmer | original | euler, k=1 | euler, k=2 | markov, k=2 | linear, k=2 |
|---|---|---|---|---|---|
| str | i64 | i64 | i64 | i64 | i64 |
| "A" | 26873 | 26873 | 26873 | 26930 | 26873 |
| "C" | 26595 | 26595 | 26595 | 26395 | 26595 |
| "G" | 19626 | 19626 | 19626 | 19564 | 19626 |
| "T" | 26906 | 26906 | 26906 | 27111 | 26906 |
| "CA" | 8444 | 7209 | 8444 | 8319 | 7815 |
| "CG" | 1321 | 5293 | 1321 | 1384 | 3369 |
| "GC" | 4843 | 5238 | 4843 | 4696 | 5050 |
| "TG" | 6735 | 5217 | 6735 | 6722 | 5997 |
All three methods keep CG rare when k=2, though in different ways: the Euler method keeps its count exactly, the Markov method approximately, and the linear method only in part, since the dinucleotides which span two blocks are not kept. Shuffling with k=1 on the other hand keeps the letter counts but brings CG back up to 5293, which is about what the letter frequencies alone would give. Only the Markov method changes the letter counts.
Why does this matter? Let us scan the original sequences and two of their shuffles for three of the example motifs, counting the hits at a P-value cutoff of 0.001. The G-box (CACGTG) and the CRE (TGACGTCA) both contain a CG, while AP-1 does not:
motifs = {m.name: m for m in um.example_motifs()}
names = ["G-box", "CRE", "AP-1"]
sets = {
"original": cg_poor,
"euler, k=1": shuffled["euler, k=1"],
"euler, k=2": shuffled["euler, k=2"],
}
pl.DataFrame({"motif": names} | {
label: [um.scan_sequences(motifs[n], s, pvalue=1e-3).height for n in names]
for label, s in sets.items()
})
| motif | original | euler, k=1 | euler, k=2 |
|---|---|---|---|
| str | i64 | i64 | i64 |
| "G-box" | 60 | 193 | 76 |
| "CRE" | 152 | 172 | 137 |
| "AP-1" | 302 | 192 | 307 |
Against the sequences shuffled with k=1, the G-box has more than three times as many hits as in the original sequences and AP-1 has a third fewer, even though nothing has changed except the order of the letters. An enrichment test against such a background would report the G-box as depleted and AP-1 as enriched, purely because of the dinucleotide composition. The background shuffled with k=2 is much closer to the original for all three motifs. For DNA, k=2 or k=3 is usually a good choice; larger values keep more of the original sequences intact, which leaves less room for shuffling, especially in short sequences.
For a single string, shuffle_string()
does the same thing. Here we can see that with k=2, the first and last
letters stay put with the Euler method:
string = "AACGTTTACGGGATCATTAG"
print(string)
print(um.shuffle_string(string, k=2, seed=1))
print(um.shuffle_string(string, k=2, seed=2))
AACGTTTACGGGATCATTAG ACGGGACATCGTTTTATAAG ACGGGAATTTTCACGTATAG
Some sequences are not uniform along their length, and a plain shuffle
would erase such local differences in composition. With window=True,
each sliding window of the sequences is shuffled in turn instead, which
keeps local patches of differing composition roughly in place. The size of
the windows and the overlap between them are set with window_size and
window_overlap, as fractions of the sequence length when below 1.
To illustrate this, let us make a sequence with an AT-rich first half and
an even second half, shuffle it both ways, and compare the AT content along
the three sequences using the windowed counts from
get_bkg():
half_at = um.create_sequences(seqnum=1, seqlen=500, seed=1,
freqs={"A": 0.4, "C": 0.1, "G": 0.1, "T": 0.4})[0]
half_even = um.create_sequences(seqnum=1, seqlen=500, seed=2)[0]
patchy = half_at + half_even
versions = {
"original": patchy,
"shuffled": um.shuffle_sequences(patchy, seed=1),
"shuffled in windows": um.shuffle_sequences(patchy, window=True, seed=1),
}
fig, ax = plt.subplots(figsize=(6, 3))
for label, s in versions.items():
at = (
um.get_bkg(s, k=1, window=True, window_size=100)
.filter(pl.col("kmer").is_in(["A", "T"]))
.group_by("start").agg(pl.col("probability").sum()).sort("start")
)
ax.plot(at["start"], at["probability"], marker="o", label=label)
ax.set_xlabel("Window start")
ax.set_ylabel("AT fraction")
ax.legend();
The plain shuffle evens out the AT content along the whole sequence, while the windowed shuffle keeps the AT-rich first half and the more even second half much as they were.
Sequence complexity¶
Low-complexity regions, such as long runs of a single letter or short
repeats, can cause trouble in motif analyses. They tend to produce
spurious motifs during de novo discovery and a pile of overlapping hits
when scanning, so it is often worth finding and removing them first.
calc_complexity() scores the complexity
of strings with one of five methods: the Wootton-Federhen and Trifonov
algorithms, a faster approximation of each, and DUST. Let us compare them
on a homopolymer, a dinucleotide repeat and a more varied string:
strings = ["AAAAAAAAAAAA", "ATATATATATAT", "ACGTTGCAAGCT"]
methods = ["WoottonFederhen", "WoottonFederhenFast", "Trifonov", "TrifonovFast", "DUST"]
pl.DataFrame({"string": strings} | {
m: um.calc_complexity(strings, m, alph="ACGT").round(3) for m in methods
})
| string | WoottonFederhen | WoottonFederhenFast | Trifonov | TrifonovFast | DUST |
|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 |
| "AAAAAAAAAAAA" | 0.0 | 0.0 | 0.0 | 0.127 | 5.0 |
| "ATATATATATAT" | 0.41 | 0.5 | 0.0 | 0.255 | 2.222 |
| "ACGTTGCAAGCT" | 0.771 | 1.0 | 0.909 | 0.982 | 0.0 |
For every method except DUST, higher values mean more complex sequences,
from 0 to 1. DUST runs the other way, with more complex sequences
approaching zero. Note that the alphabet is given here with alph, since
by default it is taken from the letters in each string, which for a
homopolymer is a single letter.
sequence_complexity() calculates
complexity in sliding windows along sequences instead, which is what we
need to find low-complexity regions. Let us take three short sequences, one
of which has a long run of As, and score them in windows of 10 with DUST:
examples = {
"A": "GTTGAAAAAAAAAAAAAAAACAGACGT",
"B": "TTAGATGGCCCATAGCTTATACGGCAA",
"C": "AATAAAATGCTTAGGAAATCGATTGCC",
}
dust = um.sequence_complexity(examples, window_size=10, method="DUST")
with pl.Config(tbl_rows=15):
display(dust)
| sequence | start | end | complexity |
|---|---|---|---|
| str | i64 | i64 | f64 |
| "A" | 0 | 10 | 0.857143 |
| "A" | 5 | 15 | 4.0 |
| "A" | 10 | 20 | 4.0 |
| "A" | 15 | 25 | 0.428571 |
| "A" | 20 | 27 | 0.0 |
| "B" | 0 | 10 | 0.0 |
| "B" | 5 | 15 | 0.0 |
| "B" | 10 | 20 | 0.0 |
| "B" | 15 | 25 | 0.0 |
| "B" | 20 | 27 | 0.0 |
| "C" | 0 | 10 | 0.285714 |
| "C" | 5 | 15 | 0.0 |
| "C" | 10 | 20 | 0.0 |
| "C" | 15 | 25 | 0.0 |
| "C" | 20 | 27 | 0.0 |
The two windows in the middle of sequence A, which cover the run of As, score 4, while every other window scores below 1. Let us keep the windows scoring at least 3:
low = dust.filter(pl.col("complexity") >= 3)
low
| sequence | start | end | complexity |
|---|---|---|---|
| str | i64 | i64 | f64 |
| "A" | 5 | 15 | 4.0 |
| "A" | 10 | 20 | 4.0 |
These regions will be masked in a moment. When whole sequences are of poor quality, it can be simpler to drop them altogether by filtering on the complexity of each sequence as a whole:
whole = um.calc_complexity(list(examples.values()), alph="ACGT")
print(whole.round(3))
{name: s for (name, s), c in zip(examples.items(), whole) if c > 0.8}.keys()
[0.601 0.862 0.817]
dict_keys(['B', 'C'])
Sequence A, with its run of As, scores lowest and is the one dropped. Where to put such a threshold is up to you, and it depends on both the method and the length of the sequences.
Sliding windows¶
Several functions above work in sliding windows, and the same machinery is
available for your own use. calc_windows()
gives the start and end of windows over a sequence of some length, counted
from 0 with the end left out as in a slice, here windows of 4 overlapping by 2 over 12 positions:
um.calc_windows(12, window=4, overlap=2)
| start | end |
|---|---|
| i64 | i64 |
| 0 | 4 |
| 2 | 6 |
| 4 | 8 |
| 6 | 10 |
| 8 | 12 |
window_string() cuts a string into those
windows, and slide_fun() applies a function
to each of them, here counting the vowels in each window:
print(um.window_string("ABCDEFGH", window=2, overlap=1))
um.slide_fun("ABCDEFGH", lambda w: sum(letter in "AEIOU" for letter in w), window=3, overlap=1)
['AB', 'BC', 'CD', 'DE', 'EF', 'FG', 'GH']
[1, 1, 1, 0]
Note that the last window can be shorter than the others; use
return_incomp=False to drop it. As a more realistic example, let us
calculate the GC content of the patchy sequence from earlier in windows of
100:
gc = um.slide_fun(patchy, lambda w: (w.count("G") + w.count("C")) / len(w), window=100)
[round(x, 2) for x in gc]
[0.23, 0.21, 0.2, 0.14, 0.2, 0.42, 0.44, 0.48, 0.49, 0.43]
Masking¶
Masking replaces parts of sequences with a placeholder letter, by default
-. Letters which are not part of the alphabet are never part of a motif
hit, so masked regions are simply skipped when scanning.
mask_seqs() masks every occurrence of a
pattern, such as runs of at least eight As:
um.mask_seqs(examples, "AAAAAAAA")
{'A': 'GTTG----------------CAGACGT',
'B': 'TTAGATGGCCCATAGCTTATACGGCAA',
'C': 'AATAAAATGCTTAGGAAATCGATTGCC'}
Overlapping occurrences are all masked, so a run of any length at least as
long as the pattern disappears entirely. rc=True masks the reverse
complement of the pattern too, and letter changes the placeholder.
mask_ranges() masks given positions
instead. These can be a list of (sequence, start, end) tuples or a
DataFrame, such as the low-complexity windows found earlier:
print(um.mask_ranges(examples, [("B", 0, 5), ("C", 19, 27)], letter="N"))
masked = um.mask_ranges(examples, low)
masked
{'A': 'GTTGAAAAAAAAAAAAAAAACAGACGT', 'B': 'NNNNNTGGCCCATAGCTTATACGGCAA', 'C': 'AATAAAATGCTTAGGAAATNNNNNNNN'}
{'A': 'GTTGA---------------CAGACGT',
'B': 'TTAGATGGCCCATAGCTTATACGGCAA',
'C': 'AATAAAATGCTTAGGAAATCGATTGCC'}
Let us check that the masked run is no longer found when scanning for a stretch of As:
poly_a = um.create_motif("AAAAAAAA", name="poly-A")
print(um.scan_sequences(poly_a, examples, pvalue=1e-3).height)
print(um.scan_sequences(poly_a, masked, pvalue=1e-3).height)
12 0
Sampling and implanting motif sites¶
A motif describes a distribution of sites, and
sample_sites() draws sites from it, choosing
each letter with the motif's probabilities at its position. For example,
let us draw a few sites from the W-box motif in the example motifs:
wbox = motifs["W-box"]
print(wbox.consensus)
um.sample_sites(wbox, n=6, seed=1)
GBTTGACYS
['TTTTGACCG', 'CGTTGACCC', 'ACTTGACTG', 'GACTGACCC', 'CTTTGTCCC', 'ACTTGACTC']
Most of the sites share the TTGAC core, while the flanks vary a fair bit
more, as they should given the motif.
implant_motifs() goes one step further
and plants sampled sites in sequences, which is useful for testing whether
an analysis can find what you know is there. With return_indices=True, it
also returns a table of where each site went:
hosts = um.create_sequences(seqnum=5, seqlen=60, seed=1)
planted, where = um.implant_motifs(wbox, hosts, seed=1, return_indices=True)
print(planted[0])
where
TTTGGACAGCGCATCTCGAATGATGTTCTTGACTAGCTTTGTCCAAAGTAGTTAAGGCGA
| sequence_i | motif_i | start | width | strand | planted |
|---|---|---|---|---|---|
| i64 | i64 | i64 | i64 | str | str |
| 0 | 0 | 26 | 9 | "+" | "TCTTGACTA" |
| 1 | 0 | 20 | 9 | "-" | "AGGTCAAAG" |
| 2 | 0 | 5 | 9 | "+" | "GCTCGACCG" |
| 3 | 0 | 37 | 9 | "+" | "GGTTGTCCG" |
| 4 | 0 | 33 | 9 | "-" | "CGGTCAAGG" |
Sites go on either strand at random, and the planted column shows each
of them as it reads on the forward strand. Several motifs can be given at
once, in which case each site comes from one of them at random. Sites can
also be placed at chosen positions, given as a list of starts, counted
from 0, for every sequence, or at a Poisson-distributed number per sequence with
rate. With centre_bias, positions are the mean of several uniform draws,
which pulls the sites towards the middle of the sequences. Let us compare
where the sites end up with and without it:
hosts = um.create_sequences(seqnum=500, seqlen=1000, seed=2)
spread = {}
for bias in (1, 4):
_, where = um.implant_motifs(wbox, hosts, centre_bias=bias, seed=1, return_indices=True)
spread[f"centre_bias={bias}"] = where["start"].quantile(0.25), where["start"].quantile(0.75)
spread
{'centre_bias=1': (276.0, 746.0), 'centre_bias=4': (381.0, 580.0)}
With centre_bias=4, the middle half of the sites starts between
positions 382 and 581, against 277 and 747 without it.
Packed sequences¶
Large sets of sequences, such as a whole genome, take up a lot of memory as
Python strings, at a byte or more per letter.
pack_sequences() packs them into a
PackedSequences object instead, where DNA
and RNA take two bits per letter:
packed = um.pack_sequences(targets)
print(packed.alphabet, len(packed), packed.nbytes)
DNA 100 25200
The 100,000 letters take up 25,200 bytes, about a quarter of a byte per
letter. Packed sequences also behave as a read-only dict, so they can go
wherever a dict of sequences can, with each sequence decoded when it is
looked up. subseq() decodes just part of a sequence:
print(packed.lengths["target_001"])
print(packed.subseq("target_001", 0, 20))
print(packed["target_001"] == targets["target_001"])
1000 TAGGGTACAATATATTTCTA True
DNA, RNA and protein letters come back in upper case, so any soft-masking
in lower case is lost. A FASTA file can be read straight into packed form
with read_fasta(path, packed=True), and
scan_sequences() scans packed sequences
without decoding them, giving the same hits as it would for the original
strings:
packed = um.read_fasta(um.example_file("fasta"), packed=True)
um.scan_sequences(wbox, packed).equals(um.scan_sequences(wbox, targets))
True
This makes packed sequences the way to go when one large set of sequences
is to be scanned many times. To scan a file just once, its path can be
passed to scan_sequences() directly
instead, in which case it is read and scanned a batch at a time.
Where next¶
With sequences and backgrounds in hand, the Scanning and enrichment guide shows how to find motif sites in them and test which motifs are enriched. The API reference has the full details of every function used here.