Genomes and signal

Genomes read from disk a piece at a time, the sequences of genomic ranges, and bigWig signal over them. Positions count from 0 and leave the end out, and ranges can be given as DataFrames or tuples, as for genomic intervals.

Genome

Genome(path: str | PathLike[str])

Bases: Mapping[str, str]

A genome in a FASTA, bgzip or 2bit file, read a piece at a time.

Plain FASTA files are read through their .fai index, and bgzipped ones through their .fai and .gzi indexes, all as written by samtools faidx; indexes that are missing are built in memory when the genome is opened, which takes a pass over the file. 2bit files carry their own index. Only the ranges asked for are read, so a genome of any size takes almost no memory until its sequences are used.

A genome is a read-only mapping from chromosome names to sequences, so it can go wherever a dict of sequences can, each chromosome being read whole when it is looked up. subseq() and extract() read just parts of them. Letters come back as they are in the file, so soft-masking is kept. scan_sequences() scans a genome one chromosome at a time, and match_bkg() reads only the windows it draws.

Examples:

>>> import universalmotif as um
>>> genome = um.Genome("TAIR10.fa.gz")
>>> genome.subseq("1", 1000, 1020)
'GGTTATAAATTATTAGAAGT'

lengths property

lengths: dict[str, int]

The length of each chromosome.

subseq

subseq(name: str, start: int, end: int) -> str

Positions start up to but not including end of the chromosome called name, counted from 0 as in a slice, reading only those.

extract

extract(ranges, nthreads: int = 1) -> dict[str, str]

The sequences of genomic ranges; see extract_ranges().

extract_ranges

extract_ranges(
    genome: Mapping[str, str], ranges, nthreads: int = 1
) -> dict[str, str]

The sequences of genomic ranges, such as peaks, from a genome.

genome is a Genome read from disk, PackedSequences or any mapping from chromosome names to sequences. ranges is a DataFrame with chrom (or seqnames, chr and the like), start and end or width columns, or (chrom, start, end[, strand]) tuples, with positions counted from 0 and the end left out, as in BED files and Python slices. Ranges on the - strand are reverse complemented. Each sequence is named chrom:start-end in the same positions, with (-) added on the minus strand, and repeated ranges get .1, .2 and so on. A genome on disk reads the ranges with nthreads threads.

Examples:

>>> import universalmotif as um
>>> genome = {"chr1": "AACCGGTTACGT"}
>>> um.extract_ranges(genome, [("chr1", 1, 5), ("chr1", 1, 5, "-")])
{'chr1:1-5': 'ACCG', 'chr1:1-5(-)': 'CGGT'}

signal_stats

signal_stats(
    bigwig: str
    | PathLike[str]
    | Mapping[str, str | PathLike[str]],
    ranges,
    stat: str = "mean",
    uncovered: float | None = 0.0,
) -> pl.DataFrame

Summarize bigWig signal, such as ATAC or ChIP coverage, over genomic ranges.

bigwig is the path of one bigWig file, giving a column called signal, or a mapping from column names to paths. ranges is a range table or (chrom, start, end[, strand]) tuples, with 0-based, half-open positions, as for extract_ranges(). stat is "mean", "max" or "min". Bases the file has no value for count as uncovered, zero by default as is usual for coverage; with uncovered=None they are left out, and a range with no values at all gets NaN.

The values are useful as covariates for match_bkg(), which can also compute them itself for the windows it draws from a genome.

Needs the bigwig extra.

Examples:

>>> import universalmotif as um
>>> peaks = [("chr1", 1000, 1200), ("chr2", 500, 700)]
>>> um.signal_stats({"atac": "atac.bw"}, peaks)

signal_values

signal_values(
    bigwig: str | PathLike[str],
    chrom: str,
    start: int,
    end: int,
    uncovered: float | None = 0.0,
) -> np.ndarray

The value of a bigWig file at each position of one genomic range.

Positions run from start up to but not including end, counted from 0. Bases the file has no value for get uncovered, zero by default, or NaN with uncovered=None, and so do all of them on a chromosome the file has no values for. Where signal_stats() summarizes signal over many ranges, this gives every position of one, such as the per-base contribution scores some models write as bigWig tracks, which can be drawn along the genome's sequence with view_attributions().

Needs the bigwig extra.

Examples:

>>> import universalmotif as um
>>> genome = um.Genome("genome.fa.gz")
>>> scores = um.signal_values("scores.bw", "1", 1000, 1100)
>>> um.view_attributions(scores, genome.subseq("1", 1000, 1100), start=1000)