# align Align deaminated sequencing reads (single- or paired-end) to a reference prepared with [`deamtools index`](index.md), writing a coordinate-sorted, indexed BAM. ## Synopsis ``` deamtools align --fasta FILE --read1 FILE --out_dir DIR --out_name NAME [options] ``` ## Required inputs | Argument | Description | |---|---| | `--fasta FILE` | Reference FASTA. Must already have been indexed with `deamtools index` (i.e. `.fai` and `.deamtools.c2t.*` exist). | | `--read1 FILE` | FASTQ for read 1, or the only FASTQ for single-end input. Plain or gzipped. | | `--out_dir DIR` | Output directory. Created automatically if it does not exist. | | `--out_name NAME` | Base name (without extension) for the output. Writes a sorted, indexed `/.bam` (and `.bam.bai`). | ## Optional arguments | Argument | Default | Description | |---|---|---| | `--read2 FILE` | *(single-end)* | FASTQ for read 2. Provide it for paired-end alignment; omit for single-end. | | `--index FILE` | *(next to the FASTA)* | Path to the converted reference built by `deamtools index` (`/.deamtools.c2t`). Use this when the index was built with a custom `--out_dir`/`--out_name`. | | `--read_group STR` | *(none)* | Read-group line passed to `bwa mem -R`, e.g. `'@RG\tID:s1\tSM:sample1\tLB:lib1\tPL:ILLUMINA'`. | | `--threads INT` | `1` | Threads used by `bwa mem` and (separately) `samtools sort`. | | `--log_level LEVEL` | `INFO` | Global flag (before the subcommand): `DEBUG`, `INFO`, `WARNING`, `ERROR`. | :::{note} By default `align` looks for the converted index next to the FASTA (`.deamtools.c2t`). If you built it with a custom location — `deamtools index --fasta ref.fa --out_dir idx --out_name foo` — point `align` at it with `--index idx/foo.deamtools.c2t`. The `index` command logs this exact path (`Converted index: …`) so you can copy it. ::: ## Requirements `bwa` and `samtools` must be on your `PATH`, and the reference must have been prepared with `deamtools index` first. ## How it works ``` FASTQ(s) ──[convert ×2]──▶ bwa mem -C ──[group + pick best]──▶ .sam 2 candidates / read | three-letter | keep higher-scoring candidate + stash original (YS) | alignment to | restore SEQ, strip f/r prefix + mark candidate (YC) | f/r conv. ref | rebuild header │ samtools sort ───┴──▶ .bam (+ .bai) ``` Every read is emitted in **two** converted forms; bwa maps both, and the better-scoring one is kept. The chosen alignments are restored and written to `/.sam`, which `samtools sort` then converts to a coordinate-sorted, indexed BAM. (The intermediate `.sam` is left in place.) ### Why two conversions per read The ACCESS-ATAC deaminase edits cytosines on **both** strands of accessible DNA. After PCR a single read can therefore show **both** `C→T` (from same-strand C deamination) **and** `G→A` (from the complementary strand's C deamination, read as a G→A relative to the forward reference). Converting a read in only one direction (the bwa-meth assumption) would leave the other direction as dense mismatches, hurting mappability for the most-edited (most accessible) reads. DeamTools instead aligns each read in **both** conversion directions and keeps whichever maps better, so a read dominated by either edit type is recovered. The reference (built by [`deamtools index`](index.md)) is doubly converted: `f` = `C→T` of the forward sequence, `r` = `G→A` of the forward sequence. ### 1. On-the-fly dual conversion As reads are streamed to the aligner, each is written twice. The two converted copies share the **original read name**, carry the original sequence in a `YS:Z:` tag, and are marked with a `YC:Z:` candidate tag (both passed through `bwa mem -C`): - **Single-end** — two candidates per read: `C→T` (`YC:Z:ct`) and `G→A` (`YC:Z:ga`). - **Paired-end** — two *fragment orientations*, with a consistent direction for the whole pair, written interleaved so `bwa mem -p` pairs each: - `f` = (read 1 `C→T`, read 2 `G→A`) - `r` = (read 1 `G→A`, read 2 `C→T`) ``` FASTQ record written to bwa: Resulting (pre-restore) SAM tags: @read1 YS:Z:ACGTTCGA YC:Z:ct SEQ = ATGTTTGA (C→T converted) ATGTTTGA YS:Z:ACGTTCGA (original, via -C) + YC:Z:ct (candidate marker) IIIIIIII ``` Reads with no quality string (e.g. FASTA input) get a placeholder quality. ### 2. Mapping and picking the best candidate `bwa mem -C -t [-p] [-R ] -` maps all candidates against the doubly-converted reference. bwa-mem preserves input order, so a read's two candidates (which share the read name) come out consecutively. Records are grouped by read name and partitioned by their `YC` tag; the candidate with the higher **primary alignment score** is kept — for pairs that is the sum of the two mates' `AS`, so a single consistent orientation (`f` or `r`) is chosen for the whole fragment. The losing candidate's records are dropped. (On a tie the first orientation, `f`/`ct`, wins.) ### 3. Restoring the alignments The kept records are rewritten back into the original reference space: - **Header.** BWA's `@SQ`/`@HD` lines describe the doubled `f`/`r` contigs, so they are dropped and replaced with a fresh `@HD VN:1.6 SO:coordinate` plus one `@SQ` per chromosome read from the original `.fai`. Other header lines (`@PG`, `@RG`, `@CO`) are passed through. - **Reference names.** The leading `f`/`r` is stripped from `RNAME` and `RNEXT`, so `fchr1`/`rchr1` both become `chr1`. - **Read sequence.** `SEQ` (the converted read) is replaced with the original from the `YS:Z:` tag — reverse-complemented for reverse-strand records (flag `0x10`), and trimmed by the CIGAR hard-clip lengths for hard-clipped records. The `YS` and `YC` tags are removed. Reference-dependent tags are restored as described below; `AS` and `XS` retain converted-reference scores. The result is a standard BAM whose coordinates, chromosome names, and read sequences are all in the original reference space, ready for `deamtools bam2bw`, `bam2fragment`, and `qc`. ### Pipeline and threading `bwa mem` runs as a subprocess: a dedicated **feeder thread** writes both converted candidates of every read into BWA's stdin while the main thread reads BWA's stdout, groups records by read name to pick the best candidate, restores them, and writes the result to `.sam`; an exception in the feeder is propagated and BWA's exit code is checked. The SAM is then converted to a coordinate-sorted BAM with `samtools sort` and indexed with `samtools index`. Both `bwa mem` and `samtools sort` use the full `--threads` count (they run one after the other). Emitting two candidates per read roughly doubles the bwa input. ## Examples ```bash # Paired-end, 8 threads deamtools align \ --fasta hg38.fa \ --read1 sample_R1.fq.gz \ --read2 sample_R2.fq.gz \ --threads 8 \ --out_dir results \ --out_name sample # -> results/sample.bam (+ results/sample.bam.bai) # Single-end with a read group deamtools align \ --fasta hg38.fa \ --read1 sample.fq.gz \ --read_group '@RG\tID:s1\tSM:sample1\tLB:lib1\tPL:ILLUMINA' \ --out_dir results \ --out_name sample # Index built in a custom location -> point align at it with --index deamtools index --fasta hg38.fa --out_dir indexes --out_name hg38 deamtools align \ --fasta hg38.fa \ --index indexes/hg38.deamtools.c2t \ --read1 sample_R1.fq.gz --read2 sample_R2.fq.gz \ --out_dir results --out_name sample ``` ## Notes - `bwa mem` and `samtools sort` each use the full `--threads` count (they run one after the other, not concurrently). - The intermediate `.sam` is kept alongside the BAM; delete it once you have the sorted BAM if you don't need it. - Paired FASTQs must contain the same number of reads in the same order; a length mismatch raises an error. - The output BAM is coordinate-sorted and indexed, so it is immediately usable by the rest of the toolkit. ### Candidate confidence and restored SAM tags Candidate selection still uses the sum of primary BWA alignment scores (`AS`) for each read or read pair. Candidates are compared in original-reference space (contig, position, strand, CIGAR and mate identity). When a competing mapped candidate has a different placement, the selected candidate's MAPQ is capped at the difference between their summed scores. Thus equal-scoring distinct placements receive MAPQ 0; duplicate conversion candidates at the same placement do not lower MAPQ. This is a conservative score-gap ceiling, not a calibrated error-probability model, and never increases BWA's MAPQ. It can reduce retention under downstream MAPQ filters compared with older DeamTools versions. `NM` and `MD` are recalculated against the original FASTA after restoring read sequences. `SA` and `XA` references and edit distances are restored as well. `AS` and `XS` remain BWA scores in converted-reference space, used for candidate selection; they are not scores recomputed against the original reference. Paired FASTQ files must have matching read names in matching order. Conventional R1 `/1` and R2 `/2` suffixes are normalized to a shared name. Mismatched names or record counts raise an error rather than producing a successful BAM. ### Per-record conversion counts DeamTools adds two custom integer tags to mapped records with restored sequence: | Tag | Definition | | --- | --- | | `ZC:i` | Number of reference C → read T observations | | `ZG:i` | Number of reference G → read A observations | Both patterns are counted regardless of the read's strand, against the original FASTA, in the same pass as NM/MD. Counts use aligned bases (`M`, `=`, `X`) without an additional base-quality filter. Insertions, deletions, reference skips and soft/hard clipping do not contribute. Zero counts are explicitly written; unmapped records or records without sequence omit both tags. Each alignment record is counted independently, including secondary/supplementary records. Mates are counted separately without overlap merging. These are observed substitutions, not proof of enzymatic editing. Do not sum these tags as a replacement for quality-filtered, mate-merged QC statistics. `NM` still includes these substitutions. `ZC`/`ZG` are DeamTools-specific conventions, not globally standardized tags; any incoming values are replaced during restoration.