Scanning large genomes¶
Scanning a few hundred promoters or peaks is quick no matter how it is done, but scanning a whole genome is a different matter. A human genome is about 3.1 billion letters long, and simply holding it in memory as Python strings takes about 3.1 GB before any scanning has even started. This guide goes through the ways the package can scan genome-sized inputs without needing that much memory: scanning straight from a FASTA file, packing sequences into a compact form which can be scanned again and again, and leaving the genome on disk to read only the parts that are needed.
A real genome would make this guide far too slow to run, so instead we will work with a small simulated one of 12 million letters. All of the timings printed below are therefore for that small example, and will be much shorter than for a real genome. At the end there are some measured figures for a human-sized genome to give an idea of what to expect at full scale.
import warnings
import polars as pl
import universalmotif as um
A small genome to work with¶
Let us make three chromosomes of three to five million letters each with
create_sequences(), with a GC content of
40% to look a little more like a real genome. Real genome assemblies also
have runs of N where the sequence is unknown, so each chromosome gets 20,000
Ns at either end, and they usually come soft-masked, meaning repeats are
written in lower case, so each chromosome also gets 40 lower case stretches
of 2,000 letters. Finally, we pick 100 random positions on each chromosome to
serve as anchors, which could be something like transcription start sites,
and plant an AP-1 site (TGACTCA) within 40 bases of most of them. These will
come in handy later on.
import numpy as np
rng = np.random.default_rng(1)
sizes = {"chr1": 5_000_000, "chr2": 4_000_000, "chr3": 3_000_000}
freqs = {"A": 0.3, "C": 0.2, "G": 0.2, "T": 0.3}
genome = {}
anchors = []
for i, (name, size) in enumerate(sizes.items()):
seq = um.create_sequences(seqnum=1, seqlen=size, freqs=freqs, seed=i + 1)[0]
for pos in np.sort(rng.integers(100_000, size - 100_000, 100)):
anchors.append((name, int(pos), int(pos) + 1))
if rng.random() < 0.6:
at = int(pos + rng.integers(-40, 41))
seq = seq[:at] + "TGACTCA" + seq[at + 7:]
seq = "N" * 20_000 + seq[20_000:-20_000] + "N" * 20_000
for start in rng.integers(20_000, size - 22_000, 40):
seq = seq[:start] + seq[start:start + 2_000].lower() + seq[start + 2_000:]
genome[name] = seq
Genomes are usually distributed as gzipped FASTA files, so let us write ours
to one in a temporary directory with
write_fasta() and Python's gzip module:
import gzip, pathlib, shutil, tempfile
tmp = pathlib.Path(tempfile.mkdtemp())
fasta = tmp / "genome.fa.gz"
um.write_fasta(genome, tmp / "genome.fa")
with open(tmp / "genome.fa", "rb") as src, gzip.open(fasta, "wb", compresslevel=1) as dst:
shutil.copyfileobj(src, dst)
Scanning straight from a FASTA file¶
The simplest way to scan a genome is to not read it in at all. When
scan_sequences() is given the path to a
FASTA file, it reads and scans the file one batch of sequence at a time, so
only a small part of the genome is ever held in memory at once. Gzipped
files, including those compressed with bgzip, are recognized from their first
bytes and decompressed as they are read. Let us scan our genome for the AP-1
motif from the example data, using two threads:
ap1 = um.example_motifs()[1]
hits = um.scan_sequences(ap1, fasta, pvalue=1e-4, nthreads=2)
hits.group_by("sequence", maintain_order=True).len()
| sequence | len |
|---|---|
| str | u32 |
| "chr1" | 1188 |
| "chr2" | 1043 |
| "chr3" | 824 |
The hits are the same as they would be for the sequences held in memory. Runs of N never form part of a hit, and lower case letters are scanned just like upper case ones, so soft-masked repeats are not skipped. Note that matches are always shown in upper case.
hits.select("sequence", "start", "end", "strand", "match", "pvalue").head()
| sequence | start | end | strand | match | pvalue |
|---|---|---|---|---|---|
| str | i64 | i64 | str | str | f64 |
| "chr1" | 24853 | 24862 | "-" | "GTGACTCAG" | 0.000038 |
| "chr1" | 27814 | 27823 | "+" | "TTGAGTCAT" | 0.000053 |
| "chr1" | 27814 | 27823 | "-" | "ATGACTCAA" | 0.000061 |
| "chr1" | 39782 | 39791 | "-" | "GTGACTCAT" | 0.000092 |
| "chr1" | 57709 | 57718 | "+" | "TTGAGTCAG" | 0.000008 |
There is one thing to be careful of here, however. The path must be given as
a path object, such as a pathlib.Path, and not as a plain string. A string
is taken to be a sequence in its own right, so scanning the text of the path
gives no hits, with a warning when the string names a file:
with warnings.catch_warnings(record=True) as caught:
print(um.scan_sequences(ap1, str(fasta), pvalue=1e-4).height)
print(caught[0].message)
0 the string '/tmp/tmpbw2eo6w8/genome.fa.gz' names a file but is scanned as a sequence; pass pathlib.Path(...) to scan the file
Packing a genome¶
Scanning from the file is ideal when a genome only needs to be scanned once.
Every scan has to read and decompress the whole file again, though, which
adds up when scanning the same genome repeatedly, for example with different
sets of motifs or P-value cutoffs. In that case it is better to read the
genome once into PackedSequences, which
store DNA and RNA in two bits per letter rather than a full byte. This can be
done straight from the file with read_fasta()
and packed=True, which packs the sequences as they are read so that the
unpacked genome is never held in memory:
packed = um.read_fasta(fasta, packed=True)
packed
PackedSequences(3 DNA sequences, 12,000,000 letters, 3.0 MB)
The packed letters take up a quarter as many bytes as the genome does as text.
print(sum(packed.lengths.values()), packed.nbytes)
12000000 3000294
Sequences which are already in memory, as a dict or a list of strings, can be
packed with pack_sequences() instead. Its
alphabet argument must match the motifs the sequences will be scanned with,
and is "DNA" by default.
um.pack_sequences(genome, alphabet="DNA", nthreads=2)
PackedSequences(3 DNA sequences, 12,000,000 letters, 3.0 MB)
What packing keeps and loses¶
Packing stores each of the four DNA or RNA letters in two bits, which leaves no room to also record their case. Upper and lower case A, C, G and T are packed the same way and come back in upper case, which means soft-masking does not survive packing. Any other letters, such as N or the IUPAC ambiguity codes, are kept separately and come back exactly as they were, case included. A short example makes this clear:
small = um.pack_sequences({"example": "ACGTNRYacgtnry"})
small["example"]
'ACGTNRYACGTnry'
Since scanning is not affected by case anyway, none of this changes the hits found in packed sequences. It only matters when the sequences are decoded again, for example to look at the soft-masking around a hit.
PackedSequences behave as a read-only
mapping from names to sequences, so they can be used anywhere a dict of
sequences can. Each sequence is only decoded when it is looked up.
print(list(packed), len(packed), "chr2" in packed)
print(packed.lengths)
print(len(packed["chr3"]))
['chr1', 'chr2', 'chr3'] 3 True
{'chr1': 5000000, 'chr2': 4000000, 'chr3': 3000000}
3000000
Looking up a whole chromosome decodes all of it, which for a real genome can
mean hundreds of megabytes at a time (and turning the whole thing into a dict
with dict(packed) undoes the savings entirely). To get a short region,
subseq() decodes only the positions
asked for, which are counted from 0 and leave the end out, as in a slice of
a string. For example, here
is the end of the run of Ns at the start of chr1, and the start of the first
soft-masked stretch, first in the original sequence and then decoded from the
packed one:
import re
print(packed.subseq("chr1", 19_990, 20_010))
masked = re.search("[acgt]", genome["chr1"]).start()
print(genome["chr1"][masked - 5:masked + 10])
print(packed.subseq("chr1", masked - 5, masked + 10))
NNNNNNNNNNGGGCAGATTG AAAGAacgcctagtt AAAGAACGCCTAGTT
Scanning packed sequences repeatedly¶
Packed sequences are scanned as they are, without decoding them first. Let us scan the genome for all twelve of the example motifs, and check that the AP-1 hits are the same as those found by scanning the file.
motifs = um.example_motifs()
all_hits = um.scan_sequences(motifs, packed, pvalue=1e-4, nthreads=2)
print(all_hits.filter(pl.col("motif") == "AP-1").drop("motif_i").equals(hits.drop("motif_i")))
with pl.Config(tbl_rows=20):
display(all_hits.group_by("motif", maintain_order=True).len())
True
| motif | len |
|---|---|
| str | u32 |
| "TATA-box" | 8071 |
| "AP-1" | 3055 |
| "CRE" | 1854 |
| "G-box" | 1349 |
| "E-box" | 1482 |
| "W-box" | 1604 |
| "GATA" | 2863 |
| "GC-box" | 470 |
| "ETS" | 1408 |
| "CArG-box" | 2869 |
| "Telo-box" | 3394 |
| "CCAAT-box" | 2142 |
The nthreads argument sets how many threads to use, and works the same way
whether the sequences come from a file, from packed sequences or from
strings. It is 1 by default, and 0 uses every core available. Let us compare
how long it takes to scan the genome for all twelve motifs with one and two
threads, both from the file and from the packed sequences:
import time
rows = []
for source, seqs in [("file", fasta), ("packed", packed)]:
for nthreads in (1, 2):
start = time.perf_counter()
um.scan_sequences(motifs, seqs, pvalue=1e-4, nthreads=nthreads)
seconds = round(time.perf_counter() - start, 2)
rows.append({"source": source, "nthreads": nthreads, "seconds": seconds})
pl.DataFrame(rows)
| source | nthreads | seconds |
|---|---|---|
| str | i64 | f64 |
| "file" | 1 | 0.69 |
| "file" | 2 | 0.4 |
| "packed" | 1 | 0.64 |
| "packed" | 2 | 0.33 |
With a genome this small every one of these is quick, and the exact times will vary from one run and one computer to the next. Scanning from the file includes the time spent decompressing and reading it, which the packed sequences only pay for once.
Packed genomes in other functions¶
Packed sequences are also useful beyond
scan_sequences(), and two functions in
particular are built to work with whole genomes.
The first is match_bkg(), which draws
background sequences matched to a set of targets on GC content and length.
With genome, it draws random windows from the genome as candidates, and
when the genome is packed only those windows are ever decoded. Let us take
200 base windows around each of our anchors as targets, as one might for a
set of peaks, and draw two matched background windows for each, keeping
clear of the targets themselves with exclude:
peaks = pl.DataFrame({
"chrom": [chrom for chrom, _, _ in anchors],
"start": [pos - 100 for _, pos, _ in anchors],
"end": [pos + 100 for _, pos, _ in anchors],
})
targets = {f"{c}:{s}-{e}": packed.subseq(c, s, e) for c, s, e in peaks.iter_rows()}
background = um.match_bkg(targets, genome=packed, exclude=peaks, n_per_target=2, seed=1)
print(len(targets), len(background))
list(background)[:3]
300 600
['chr1:1557449-1557649', 'chr3:63506-63706', 'chr1:357203-357403']
The targets and background can then be compared, for example with
enrich_motifs().
enriched = um.enrich_motifs(motifs, targets, background, nthreads=2)
enriched.select("motif", "consensus", "target_seq_hits", "bkg_seq_hits", "enrichment", "qvalue")
| motif | consensus | target_seq_hits | bkg_seq_hits | enrichment | qvalue |
|---|---|---|---|---|---|
| str | str | i64 | i64 | f64 | f64 |
| "AP-1" | "NTGASTCAN" | 182 | 15 | 24.266667 | 2.1431e-87 |
| "W-box" | "GBTTGACYS" | 24 | 20 | 2.4 | 0.013033 |
The planted AP-1 sites stand out clearly, as they should. The W-box motif is also enriched, though to a much lesser degree, since its TTGAC core overlaps the planted TGACTCA whenever the base before it happens to be a T.
The second is motif_proximity(), which
tests whether motif hits sit closer to a set of anchors than expected by
chance. Given motifs and a packed genome, it scans the genome itself without
decoding it.
proximity = um.motif_proximity(anchors=anchors, motifs=motifs, sequences=packed, nthreads=2)
proximity.select("motif", "nhits", "best_dist", "hits_near", "expected_near", "enrichment", "qvalue")
| motif | nhits | best_dist | hits_near | expected_near | enrichment | qvalue |
|---|---|---|---|---|---|---|
| str | i64 | i64 | i64 | f64 | f64 | f64 |
| "AP-1" | 1704 | 50 | 184 | 4.3026 | 42.76484 | 2.0101e-226 |
| "W-box" | 1604 | 50 | 21 | 4.0501 | 5.185057 | 1.5739e-7 |
Again the AP-1 hits are found to be concentrated close to the anchors, with
the most significant enrichment within 50 bases of them, which fits the sites
having been planted within 40 bases. The W-box motif shows up for the same reason as
before. The default binomial test assumes hits are spread evenly across the
genome, so for real data it is worth also trying method="permutation",
which is slower but accounts for biases in sequence composition.
Leaving the genome on disk¶
Packing still means reading the whole genome once, and holding a quarter of
it in memory afterwards. Often only small parts of a genome are needed, such
as the sequences of a set of peaks or a few thousand background windows, and
for those a Genome is the better fit. It leaves
the genome where it is and reads only the parts asked for, using an index of
where each sequence starts in the file.
Genomes can be opened from any of three formats: plain FASTA files, FASTA
files compressed with bgzip, and 2bit files. A bgzipped file is compressed
in small blocks that can each be decompressed on their own, which is what
lets a piece of it be read without decompressing everything before it. The
one format that cannot be read this way is a FASTA file compressed with
ordinary gzip, which is how genomes are often downloaded, and which is also
how our genome.fa.gz was written. Opening such a file gives an error
saying so, and the fix is to recompress it with bgzip, which comes with
samtools and htslib:
gunzip genome.fa.gz
bgzip genome.fa # writes genome.fa.gz again, now in blocks
samtools faidx genome.fa.gz
The last step writes the .fai and .gzi indexes next to the file. They
are optional, but without them the indexes are built each time the genome is
opened, which means reading through the whole file once. With them, opening
a genome takes no time at all. The same goes for the .fai index of a plain
FASTA file, while 2bit files carry their own index.
Since we have no bgzip here, let us use the uncompressed genome.fa
written along the way. A bgzipped or 2bit genome is opened in exactly the
same way, with the format recognized from the file itself.
on_disk = um.Genome(tmp / "genome.fa")
on_disk
Genome('genome.fa', 3 sequences, 12,000,000 letters)
A genome on disk is a read-only mapping from names to sequences, just like
packed sequences, and subseq() reads a
region counted from 0 and leaving the end out. Letters come back exactly as
they are in the file, so the soft-masking that packing loses is kept:
print(on_disk.lengths)
print(on_disk.subseq("chr1", masked - 5, masked + 10))
{'chr1': 5000000, 'chr2': 4000000, 'chr3': 3000000}
AAAGAacgcctagtt
The sequences of a whole set of ranges can be read at once with
extract_ranges(), which takes the same
range tables as match_bkg() and
motif_proximity(). Ranges on the minus
strand come back reverse complemented. Here are the sequences of our peaks,
which match the targets taken from the packed genome once the case is set
aside:
peak_seqs = um.extract_ranges(on_disk, peaks, nthreads=2)
print(list(peak_seqs) == list(targets))
print(all(a.upper() == b for a, b in zip(peak_seqs.values(), targets.values())))
um.extract_ranges(on_disk, [("chr1", masked - 5, masked + 10, "-")])
True True
{'chr1:103424-103439(-)': 'aactaggcgtTCTTT'}
Hits found by scanning these sequences are counted from the start of each
peak, which is rarely what is wanted afterwards.
genomic_hits() puts them back on the
genome, reading the position of each peak from its name, and flips the
strand of any hit in a peak that was reverse complemented:
peak_hits = um.scan_sequences(ap1, peak_seqs, pvalue=1e-4)
on_genome = um.genomic_hits(peak_hits)
on_genome.select("sequence", "chrom", "start", "end", "strand", "match").head()
| sequence | chrom | start | end | strand | match |
|---|---|---|---|---|---|
| str | str | i64 | i64 | str | str |
| "chr1:267192-267392" | "chr1" | 267295 | 267304 | "+" | "ATGACTCAG" |
| "chr1:267192-267392" | "chr1" | 267295 | 267304 | "-" | "CTGAGTCAT" |
| "chr1:289945-290145" | "chr1" | 290068 | 290077 | "+" | "CTGACTCAG" |
| "chr1:289945-290145" | "chr1" | 290068 | 290077 | "-" | "CTGAGTCAG" |
| "chr1:359726-359926" | "chr1" | 359802 | 359811 | "+" | "TTGACTCAG" |
Every one of these is among the hits found earlier by scanning the whole genome from the file:
key = ["chrom", "start", "end", "strand"]
print(on_genome.height, on_genome.join(hits.rename({"sequence": "chrom"}), on=key, how="semi").height)
330 330
Ranges can also come from, and go back to, the genomic interval libraries.
PyRanges objects and polars-bio frames can be given anywhere ranges are
taken, since they count positions the same way, or are converted when a
polars-bio frame is marked as counting from 1, and
to_pyranges() and
to_polars_bio() hand tables such as these
hits back to them, for example to find the hits overlapping a set of genes:
import polars_bio as pb
genes = pb.read_gff("genes.gff3")
pb.overlap(um.to_polars_bio(on_genome), genes)
match_bkg() reads only the windows it draws
from a genome on disk, and with the same seed it draws the same windows as
from the packed genome. motif_proximity()
reads and scans one chromosome at a time, and so does
scan_sequences() when given a genome on
disk, so only a single chromosome is ever held in memory.
background_disk = um.match_bkg(targets, genome=on_disk, exclude=peaks, n_per_target=2, seed=1)
print(list(background_disk) == list(background))
proximity_disk = um.motif_proximity(anchors=anchors, motifs=motifs, sequences=on_disk, nthreads=2)
proximity_disk.drop("dist_to_nearest").equals(proximity.drop("dist_to_nearest"))
True
True
Matching backgrounds on signal¶
GC content and length are not the only things that set peaks apart from the
rest of the genome. Peaks from ChIP-seq or ATAC-seq also sit in open
chromatin, and open chromatin has its own sequence character, so a
background drawn from anywhere in the genome can make motifs look enriched
simply because they are common in open regions. If the signal is available
as a bigWig file, match_bkg() can match on it
as well: give the file with bigwig and the coordinates of the targets
with target_ranges, and the mean signal over each target and each window
becomes one more value to match on.
background = um.match_bkg(
targets, genome=on_disk, exclude=peaks, seed=1,
bigwig={"atac": "atac.bw"}, target_ranges=peaks,
)
The signal over any set of ranges can also be had directly with
signal_stats(), for example to check how
well the background matches. Both need the bigwig extra, installed with
pip install universalmotif[bigwig].
Memory and speed in practice¶
To give an idea of what to expect from a real genome, here are some figures measured on a human-sized genome of 3.1 billion letters, using two threads. These are approximate, and will depend on the computer, the motifs and the P-value cutoff:
- Scanning one motif straight from the FASTA file took about 5.5 seconds and about 0.25 GB of memory. Reading a gzipped file added a few seconds for decompression.
- Scanning 20 motifs straight from the file at a P-value cutoff of 1e-5 took about 70 seconds.
- Packing the genome took about 3.3 seconds, and the packed genome held about 0.8 GB, compared to about 3.1 GB for the same genome held as Python strings.
- Opening an indexed genome on disk took no time at all, and reading 2,000 random regions of up to 3,000 letters from a bgzipped plant genome took about 0.1 seconds.
So, in general, if a genome only needs to be scanned once, scanning it
straight from the file is the way to go, as memory stays low regardless of
the size of the genome. If it is going to be scanned several times, or used
with match_bkg() or
motif_proximity(), packing it once is
usually worth the extra memory. When only parts of the genome are needed,
such as the sequences of peaks or background windows, a
Genome on disk needs barely any memory at all.
Reading a genome into a dict of strings is rarely necessary.
Note that the hits themselves also take memory. By chance alone, a cutoff of
1e-4 gives about one hit for every 10,000 positions on each strand, which for
a 3.1 billion letter genome works out to around 600,000 hits per motif. For
genome-wide scans a stricter cutoff such as 1e-5 or 1e-6 is usually more
sensible, and no_overlaps=True can further remove redundant overlapping
hits.
A stricter cutoff comes with a caveat, however. A short motif only has so
many possible matches, and its very best match may simply not be rare enough
to reach a strict cutoff. The lowest P-value a motif can reach is that of its
maximum score, which we can get with
motif_range() and
motif_pvalue():
lowest = pl.DataFrame({
"motif": [m.name for m in motifs],
"width": [m.width for m in motifs],
"lowest_pvalue": [um.motif_pvalue(m, score=um.motif_range(m)[1]) for m in motifs],
})
with pl.Config(tbl_rows=20):
display(lowest)
| motif | width | lowest_pvalue |
|---|---|---|
| str | i64 | f64 |
| "TATA-box" | 10 | 9.5367e-7 |
| "AP-1" | 9 | 0.000008 |
| "CRE" | 8 | 0.000015 |
| "G-box" | 8 | 0.000015 |
| "E-box" | 8 | 0.000015 |
| "W-box" | 9 | 0.000004 |
| "GATA" | 8 | 0.000031 |
| "GC-box" | 9 | 0.000004 |
| "ETS" | 10 | 9.5367e-7 |
| "CArG-box" | 10 | 9.5367e-7 |
| "Telo-box" | 9 | 0.000004 |
| "CCAAT-box" | 9 | 0.000004 |
None of the four eight letter motifs can reach a P-value below 1e-5 at all, so scanning with that cutoff will never give them a single hit, no matter how large the genome. Likewise, a cutoff of 1e-6 leaves only the ten letter motifs able to produce any hits. It is worth checking this before settling on a cutoff for a set of motifs of varying widths.
Finally, on the matter of threads, more threads will generally make scanning faster, though on a shared computer it is polite to leave some cores free for others. The figures above were measured with only two threads, which is usually plenty to scan a genome in reasonable time.
Where next¶
The Gapped motifs and other alphabets guide covers scanning with gapped, RNA, amino acid and custom alphabet motifs, and the API reference has the details of every function used here.