# 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`

```python
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`

```python
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`

```python
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`

```python
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`

```python
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](https://github.com/ENCODE-DCC/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`

```python
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`

```python
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`

```python
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`

```python
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 |
|---|---|
| `get_chrom_sizes_from_bam(bam)` | Chromosome name → length, from an open BAM header. |
| `get_chrom_sizes_from_file(path)` | Chromosome name → length, from a UCSC `chrom.sizes` file. |
| `get_version()` | Installed package version (via `importlib.metadata`). |
| `_load_regions(bed_path)` | BED as a `chrom/start/end` DataFrame, with overlapping intervals merged. |