# Algorithm This page describes how DeamTools aligns deaminated reads and detects, quantifies, and summarises deamination events. ## Biological basis Double-stranded DNA cytosine deaminases convert cytosine (C) to uracil (U) on DNA that is not shielded by nucleosomes or bound proteins. Accessible chromatin is therefore edited at a high rate, while protein-bound DNA is protected, leaving **footprints**. After PCR (which reads U as T) and sequencing, a deaminated cytosine appears as a C→T substitution relative to the reference; on the opposite strand the same event appears as G→A. The signal is a mismatch: the reference carries C (or G), but the aligned read carries T (or A). The density of these mismatches reports accessibility, and local depletions reveal bound factors. ## Strand convention The two mismatch patterns say **which strand of the molecule was deaminated** — not which way the observing read happened to align. **Top-strand event.** A C on the plus strand is deaminated: ``` Reference (+): 5'—...C...—3' ↓ deamination Molecule (+): 5'—...T...—3' ``` Pattern: `ref_base == 'C'` and `read_base == 'T'`. **Bottom-strand event.** A C on the minus strand is deaminated. That C pairs with a reference G, so in reference coordinates the event shows up at a G: ``` Reference (+): 5'—...G...—3' Template (−): 3'—...C...—5' ← deaminated here ``` Pattern: `ref_base == 'G'` and `read_base == 'A'`. **Read orientation does not select between them.** Two facts make `is_reverse` the wrong test: 1. A BAM stores `SEQ` in **reference orientation** — a reverse-strand record has already been reverse-complemented — so both patterns are expressed against the plus strand regardless of how the read aligned. 2. Library prep fixes the deaminated U into a real **T:A** base pair, so after amplification *both* strands of the molecule carry the substitution and any read covering the position reports it. A double-stranded deaminase edits both strands, so one fragment routinely carries both patterns at different positions. Every command therefore calls both patterns on every read. Measured on `data/ACCESS-ATAC/chr10.bam`, `is_reverse` predicts a read's dominant edit direction only about **65%** of the time; filtering on it discarded **35.5%** of editing events — events carrying the same `TC` enzyme-motif fingerprint as the ones kept, so genuine deamination rather than noise. `bam2fragment` used such a filter until 2026-09-10. ## Deamination-aware alignment (`index` + `align`) Heavily edited reads align poorly to an unmodified reference, so DeamTools uses a bwa-meth-style three-letter strategy. - **`index`** writes a *doubly-converted* reference: every chromosome appears twice, once with all C→T (prefixed `f`) and once with all G→A (prefixed `r`), then runs `bwa index` on it. - **`align`** maps each read in **both** conversion directions and keeps whichever scores higher. Because the deaminase edits both strands, a single read can carry both `C→T` and `G→A` edits, so each read is emitted twice: single-end as `C→T` (`YC:Z:ct`) and `G→A` (`YC:Z:ga`); paired-end as two fragment orientations `f` = (R1 `C→T`, R2 `G→A`) and `r` = (R1 `G→A`, R2 `C→T`). The original sequence is stashed in `YS:Z:` and the candidate in `YC:Z:`, both carried through `bwa mem -C`. After mapping, records are grouped by read name and the candidate with the higher primary alignment score (sum of the mates' `AS` for pairs) is kept; the original SEQ is restored from `YS`, the `f`/`r` prefix is stripped from RNAME/RNEXT, the `YS`/`YC` tags are dropped, and the BAM is sorted and indexed. Choosing one orientation per fragment keeps the mates on the same converted contig, so proper pairing is preserved. ### Read space vs reference space `f` = (R1 `C→T`, R2 `G→A`) reads as though the two mates were converted differently. **They are not.** Each candidate is *one* conversion applied to the whole fragment in **reference space**; the labels differ only because R2's raw FASTQ sequence is the reverse complement of the reference-space sequence, and the conversion has to be written in the space the FASTQ is actually in. A fragment whose **top** strand was deaminated, with R1 reading its left end and R2 its right end: ``` reference (+) ACGTCGATCGTTAGCCATGC deaminated at * * * molecule, top ATGTCGATTGTTAGCTATGC R1 (read space) ATGTCGATTG ← already reference space R2 (read space) GCATAGCTAA ← reverse complement of it R2 in ref space TTAGCTATGC ``` Candidate `f` means "convert this fragment `C→T` in reference space": ``` R1 read-space C→T ATGTTGATTG R2 read-space G→A ACATAACTAA ──bwa reverse-complements──▶ TTAGTTATGT R2 ref-space C→T TTAGTTATGT ✓ identical ``` The two agree because complementing turns `G→A` into `C→T`. Formally, for any sequence `s`: ``` revcomp(s.replace("G", "A")) == revcomp(s).replace("C", "T") ``` So `f` puts **both** mates on the `f` contig and `r` puts both on `r`. That is the whole point of the label flip: it is what keeps the pair on one contig. ### Why there are two candidates and not four Converting each mate independently would give four combinations. The two extra ones are not alternative hypotheses about the molecule — they place the mates on *different* contigs, because a `C→T` mate goes to `f` and a `G→A` mate to `r`. Measured on 300 simulated top-strand-deaminated fragments against `data/test/chr20sim.deamtools.c2t`: | candidate | mates on one contig | proper pairs | mean `AS(R1)+AS(R2)` | |---|---|---|---| | `f` = (R1 `ct`, R2 `ga`) | 278/279 | 273/279 | **200.0** | | `r` = (R1 `ga`, R2 `ct`) | 272/279 | 226/279 | 159.2 | | (R1 `ct`, R2 `ct`) | 20/279 | **0/279** | 177.4 | | (R1 `ga`, R2 `ga`) | 19/279 | **0/279** | 176.9 | Each mate still aligns somewhere on its own, so the invalid combinations collect a respectable `AS` — 177.4 here, **above** the valid `r` candidate's 159.2. Adding them to the take-best would therefore sometimes select a cross-contig, non-proper pair over a correct one, and after the `f`/`r` prefix is stripped the two mates would land on the same chromosome name carrying mate coordinates and `TLEN` that mean nothing. The reason only one choice exists per fragment is biological, not a simplification: R1 and R2 report the *same molecule* in the same reference orientation, so whatever edits it carries appear in both. The conversion is a property of the fragment, so there is one choice to make, and take-best makes it. What take-best cannot do is absorb **both** directions at once — a double-stranded deaminase edits both strands, so a real fragment carries `C→T` and `G→A` together and neither three-letter alphabet covers both. That is a genuine limit of the three-letter approach, and it is why `deamtools` accuracy declines gently with editing rate (99.26% → 99.11% from 0% to 100% editing) rather than staying flat. Fixing *that* would need a different alignment model, not more conversion candidates. ## Edit detection conventions All commands iterate aligned positions with `pysam`'s `get_aligned_pairs(matches_only=True)`, so only matched (M/=/X) bases are compared — insertions, deletions, and clips are skipped, and query/reference bases are always directly comparable. Two strand conventions are used, depending on the command: | Command | Convention | Counted as an edit | |---|---|---| | `bam2bw`, `qc` | **strand-agnostic** | any reference `C→T` **or** `G→A` mismatch, regardless of read orientation | | `bam2fragment` | **strand-agnostic, strand-resolved** | the same mismatches, reported in two separate columns so the edited strand is preserved | Reads flagged unmapped, duplicate, QC-fail, secondary, or supplementary are always excluded, then `min_mapq` is applied per read; `min_baseq` gates individual bases. ### Mates are merged before counting `bam2bw` and `qc` count per **fragment**, not per alignment record. Whenever the insert is shorter than twice the read length the two mates overlap, and every reference position in that overlap is reported by both of them — one position on one molecule, seen twice. Records are therefore grouped by read name and each fragment's mates are collapsed into one `reference position → base` map before anything is tallied, so an overlapping position contributes once. Where the mates disagree at an overlap, the **higher base quality wins**: library prep turns the deaminated C into a real T:A pair that both strands carry, so the two mates are expected to agree and a disagreement means one of them misread. This is not a rounding detail. On `data/ACCESS-ATAC/chr10.bam` the overlap accounted for **24.5%** of `qc`'s editing opportunities and **21.2%** of `bam2bw`'s count-mode signal. Because the overlap is the middle of the fragment rather than a random subset of positions, counting it twice biased every rate derived from it rather than merely inflating the totals. Records whose mate never arrives — unpaired reads, a mate that was unmapped, filtered out, or placed on another contig — are counted as one-record fragments, so nothing is dropped. Single-end data is unaffected: one record is one fragment. `bam2fragment` has always worked per fragment, and now shares the same merging helper. ## Tn5 cut sites (`bam2bw --event tn5`) Deaminase libraries such as ACCESS-ATAC are tagmented with Tn5, so each fragment end also marks a Tn5 insertion — the ATAC-style accessibility signal. With `--event tn5`, `bam2bw` counts those insertion sites instead of edits. Each passing read contributes the insertion at its **5′ end**. For paired-end data the two mates' 5′ ends are the two ends of the fragment, so every fragment gives two cuts; for single-end data each read gives the one cut it was sequenced from. A reverse read's 5′ end is the right end of its alignment, not its left-most coordinate. Tn5 inserts as a dimer that nicks the two strands 9 bp apart, so the fragments on either side of one insertion share a 9-bp duplication `[p, p+9)`: ``` reference ......[p ....... p+8]...... right fragment, fwd [p ─────────────────▶ 5′ base = p left fragment, rev ◀────────────── p+8] 5′ base = p+8 (reference_end = p+9) shifted cut p+4 ``` Moving each read 4 bp inward from its 5′ base — `reference_start + 4` forward, `reference_end − 5` reverse, the usual **+4/−5 shift** — puts both reads of one insertion on the same base, `p+4`, the centre of the duplication. Both shifts are options — `--forward_shift` (default `+4`, added to the start) and `--reverse_shift` (default `−5`, added to the exclusive end) — so a different convention, a pre-shifted BAM (`0`/`0`), or raw 5′ ends (`0`/`−1`) are one flag away. No reference sequence is involved, so `--fasta` is optional in this mode. `min_baseq` does not apply (a cut has no base to be of poor quality); the flag filters, `min_mapq`, `extend_size` and `normalize` do. `--mode ratio` is edit-only and is rejected with `--event tn5`. The geometry can be checked on real data by cross-correlating forward and reverse 5′ bases: pairs from one insertion should sit at an offset of +8. On the Ultima single-end concurrent ACCESS-ATAC BAM they do, strongly (5.7× background at +8 against 2.4× at +9). On the Illumina paired-end `HepG2.bam` there is a second population of similar size at +9 (2.3× at +8, 2.6× at +9), for which the shifted cuts land 1 bp apart. It is not soft-clipping (excluding clipped 5′ ends changes nothing) and not an extra non-genomic 5′ base (the first base's mismatch rate is 0.28%, no higher than the next); its origin is not yet known. ## Signal generation (`bam2bw`) For each region (whole chromosome or a merged BED interval) the reference is fetched once and reads are streamed via the BAM index. - **count mode** — a per-base count of edits. With `--extend_size E > 0`, each event is broadcast symmetrically into a window of width `2E + 1` (clipped to the region), so the value at a base is the number of events within `E` bp. - **ratio mode** — `edit_count / total_ACGT_coverage` at each base; positions whose total coverage is below `--min_coverage` are written as `0`. `--extend_size` is ignored in this mode. Both are per fragment. The denominator is counted in the same pass as the numerator, over the same merged fragments, so the two cannot disagree about which reads they saw — and it now honours `--min_mapq`. Before 2026-09-10 the denominator came from pysam's `count_coverage`, which applies no mapping-quality filter, so ratio mode was dividing edits from MAPQ-passing reads by coverage that included reads the numerator had excluded. BED intervals are merged before counting so a read spanning an overlap is not double-counted: ``` Input: [0, 500) [300, 800) [1000, 1500) Merged: [0, 800) [1000, 1500) ``` ## Enzyme sequence bias Deaminases have intrinsic flanking-sequence preferences (e.g. DddA strongly prefers `TC`). `qc` reports the per-trinucleotide edit fraction and renders a deaminase sequence-motif logo (the reference window around each editing event, centre excluded, `G→A` events reverse-complemented), both built directly from the BAM. A strongly skewed profile indicates that footprinting/occupancy analyses should correct for enzyme bias. ## Parallelism and output `bam2bw` first groups its regions into **batches** of neighbouring intervals (about eight per worker, capped at 50 Mbp of reference each), and each batch opens the BAM once. Opening a BAM reads its whole index, so one open per region used to dominate peak runs: 57k HepG2 peaks meant 1.6–3.6 minutes of index loading alone. Each chromosome (or batch) is processed in its own **worker process** (`concurrent.futures.ProcessPoolExecutor`, via `deamtools.utils.run_jobs`); BAM and FASTA handles are opened independently per worker. Processes rather than threads because the per-read loops are pure Python, so a thread pool is serialised by the GIL — `qc --threads 12` on a 3.2 M-record BAM ran at 105% CPU as threads and runs about 6× faster as processes. Results are merged in submission order, not completion order, so the output is bit-for-bit identical whatever the worker count. With `--threads 1` everything runs in the calling process. For BigWig output, regions are sorted into `(chrom, start, end)` order — as required by `pyBigWig` — and only non-zero bases are written (sparse, `span=1`); the header lists every chromosome even if it has no signal. ## Complexity For *N* aligned reads of average length *L*, edit detection is O(*N* × *L*) per thread. `bam2bw` holds one float32 array per in-flight chromosome (≈ one float per base), so peak memory scales with the largest chromosomes times the number of threads.