API reference
Every subcommand is a thin wrapper around a public run_* function, so DeamTools can be used as a Python library as well as from the command line. Each function is importable from its submodule.
Note
With threads > 1, run_bam2bw, run_bam2fragment, run_qc and run_footprint run their work in worker processes. On macOS, where Python starts workers by re-importing your script, keep the call under the usual guard:
if __name__ == "__main__":
run_qc(..., threads=8)
Without it each worker re-runs the script, and Python stops with a RuntimeError about the bootstrapping phase. Jupyter notebooks need no guard. Code piped in on stdin (python -) cannot be re-imported at all, so the work runs in a single process there, with a warning.
deamtools.align.index
run_index
from deamtools.align.index import run_index
run_index(
fasta_path: str,
out_dir: str | None = None,
out_name: str | None = None,
force: bool = False,
) -> None
Build the deamination-aware reference index: the standard <fasta>.fai (always next to the FASTA) plus the doubly-converted reference and its BWA index at <out_dir>/<out_name>.deamtools.c2t*. out_dir/out_name default to the FASTA’s own directory and file name. Requires bwa and samtools.
deamtools.align.align
run_align
from deamtools.align.align import run_align
run_align(
fasta_path: str,
read1: str,
out_dir: str,
out_name: str,
read2: str | None = None,
threads: int = 1,
read_group: str | None = None,
index_path: str | None = None,
) -> None
Align deaminated reads (read 1 C→T, read 2 G→A) to a deamtools index-prepared reference and write a sorted, indexed <out_dir>/<out_name>.bam. Pass index_path when the index was built with a custom --out_dir/--out_name. Requires bwa and samtools.
deamtools.preprocessing.bam2bw
run_bam2bw
from deamtools.preprocessing.bam2bw import run_bam2bw
run_bam2bw(
bam_path: str,
fasta_path: str | None, # required for event="edit"
out_dir: str,
out_name: str,
chrom_sizes_path: str | None = None,
bed_path: str | None = None,
min_mapq: int = 20,
min_baseq: int = 20,
extend_size: int = 0,
threads: int = 1,
mode: str = "count",
min_coverage: int = 10,
normalize: bool = False,
scale_factor: float = 1_000_000,
event: str = "edit", # or "tn5" for Tn5 insertion sites
forward_shift: int | None = None, # tn5: default +4
reverse_shift: int | None = None, # tn5: default -5 (on the exclusive end)
) -> None
Write a per-base BigWig of deamination signal to <out_dir>/<out_name>.bw. mode="count" writes raw edit counts (strand-agnostic C→T or G→A) and honours extend_size; mode="ratio" writes edits / total_ACGT_coverage, masking positions below min_coverage to 0. In count mode, normalize=True scales every value by scale_factor / total so the track sums to scale_factor (reads/counts-per-million); ignored in ratio mode.
Both numerator and denominator are counted per fragment: the mates of a pair are merged first, so a reference position both mates cover contributes once rather than twice, and the higher base quality settles a disagreement. They are computed in the same pass over the same fragments and both honour min_mapq.
fasta_path may be None with event="tn5". event="tn5" counts Tn5 insertion sites instead: each passing read’s 5′ end, shifted +4 (forward) / −5 from the exclusive end (reverse), so a paired-end fragment contributes both ends and a single-end read its start. Count mode only.
deamtools.preprocessing.bam2fragment
run_bam2fragment
from deamtools.preprocessing.bam2fragment import run_bam2fragment
run_bam2fragment(
bam_path: str,
fasta_path: str,
out_dir: str,
out_name: str,
min_mapq: int = 20,
min_baseq: int = 20,
threads: int = 1,
barcode: bool = False,
barcode_tag: str = "CB",
gzip: bool = False,
) -> None
Write a per-fragment editing-signal table to <out_dir>/<out_name>.tsv (or .tsv.gz when gzip=True). With barcode=True, a barcode column is added (10x ordering). The last two columns are the C→T and G→A positions, reported separately: C→T means the top strand was deaminated at that position and G→A the bottom strand. Both are called regardless of read orientation, and a position covered by both mates is reported once.
deamtools.qc
run_qc
from deamtools.qc import run_qc
run_qc(
bam_path: str,
fasta_path: str,
out_dir: str,
out_name: str,
tss_path: str | None = None,
min_mapq: int = 20,
min_baseq: int = 20,
threads: int = 1,
tss_flank: int = 2000,
plot: bool = True,
n_reads: int | None = None,
logo_scale: str = "frequency", # or "bits"
) -> dict
Compute QC metrics and write <out_dir>/<out_name>.json plus a self-contained HTML report. The report includes the deaminase sequence-motif logo, built directly from the editing events in the BAM. Returns the metrics dictionary.
Editing is counted per fragment: for paired-end data the mates are merged first, so a reference position both mates cover is one observation rather than two, and the higher base quality decides a disagreement. The reads block still counts records; fragments reports how many fragments the editing metrics ran over and how many came from merged pairs.
n_reads subsamples the BAM to roughly that many reads (drawn uniformly across the genome from a fixed seed) instead of reading all of them; rates stay unbiased, but the reported counts are counts of the sample. The draw is keyed on the read name, so both mates of a fragment always share its fate.
Supplying tss_path adds TSS enrichment, computed the way the ENCODE ATAC-seq pipeline defines it: insertion 5’ ends binned at 10 bp over a ±tss_flank window, minus-strand TSS flipped, normalised to the outermost 100 bp on each side, and scored as the peak of that profile. The report gains a TSS plot with the score marked on it, and the per-bin numbers behind it are written to <out_dir>/<out_name>.tss_enrichment.csv.
deamtools.motif.match
run_motif_matching
from deamtools.motif.match import run_motif_matching
run_motif_matching(
fasta_path: str,
bed_path: str,
out_dir: str,
out_name: str,
motifs: list | None = None,
release: str = "JASPAR2024",
collection: str = "CORE",
tax_group: list[str] | None = None,
pseudocounts: float = 0.0001,
p_value: float = 1e-4,
) -> None
Scan the sequence of each BED region with motifmatchpy and write motif matches to <out_dir>/<out_name>.bed as 6-column BED (chrom, start, end, motif, score, strand). Motifs are fetched from JASPAR (needs pyjaspar) unless motifs is passed explicitly.
prepare_scanner / scan_sequence
from deamtools.motif.match import prepare_scanner, scan_sequence
scanner = prepare_scanner(motifs, pseudocounts=0.0001, p_value=5e-05)
matches = scan_sequence(scanner, motifs, seq, chrom, offset=0)
load_motifs_from_files reads .pfm/.adm files into motifmatchpy.Motif objects, naming each after its file stem. prepare_scanner builds a motifmatchpy.MotifScanner (log-odds matrices, p-value thresholds, both strands). scan_sequence scans one sequence and returns (chrom, start, end, name, score, strand) tuples. These let you scan in-memory sequences/motifs without writing a BED.
deamtools.footprint
run_footprint
from deamtools.footprint import run_footprint
run_footprint(
bigwig_path: str,
regions_path: str,
out_dir: str,
out_name: str,
n_shuffles: int = 1000,
threads: int = 1,
seed: int | None = None,
) -> None
Score TF footprints at motif sites. For each site of width L, reads the per-base BigWig over the 3 * L window [start - L, end + L) and computes fp_score = mean(left flank) + mean(right flank) - mean(centre); positive scores get a permutation p-value (n_shuffles within-window permutations). Writes <out_dir>/<out_name>.bed with columns chrom, start, end, name, fp_score, p_value. Sites whose window falls off the chromosome are skipped.
deamtools.utils
from deamtools.utils import (
get_chrom_sizes_from_bam, # (bam: pysam.AlignmentFile) -> dict[str, int]
get_chrom_sizes_from_file, # (path: str) -> dict[str, int]
get_version, # () -> str
)
from deamtools.utils.regions import _load_regions # (bed_path: str) -> pandas.DataFrame
Function |
Returns |
|---|---|
|
Chromosome name → length, from an open BAM header. |
|
Chromosome name → length, from a UCSC |
|
Installed package version (via |
|
BED as a |