# bam2bw

Convert a coordinate-sorted BAM file to a BigWig track of per-base deamination counts — or, with `--event tn5`, of Tn5 insertion (cut) sites.

Signal is counted **per fragment**, not per alignment record: for paired-end data the two mates are merged before counting, so a reference position that both mates cover contributes once instead of twice, and the higher base quality settles a disagreement between them. In `--mode ratio` the denominator is counted the same way, in the same pass, and honours `--min_mapq`. See [Algorithm](../algorithm.md#mates-are-merged-before-counting) for why this matters — on real ACCESS-ATAC data the mate overlap was 21% of the count-mode signal.

## Synopsis

```
deamtools bam2bw --bam FILE --fasta FILE --out_dir DIR --out_name NAME [options]
```

## Required inputs

| Argument | Description |
|---|---|
| `--bam FILE` | Coordinate-sorted BAM file. Must be accompanied by an index (`.bai`). |
| `--fasta FILE` | Reference FASTA file used during alignment. Must be indexed with `samtools faidx` (`.fai`). **Required with `--event edit`** (the default). With `--event tn5` it is optional and never opened, since a cut site depends only on where a read aligns. |
| `--out_dir DIR` | Output directory. Created automatically if it does not exist. |
| `--out_name NAME` | Base name (without extension) for the output. The BigWig is written to `<out_dir>/<out_name>.bw`. |

## Optional arguments

### Input / scope

| Argument | Default | Description |
|---|---|---|
| `--chrom_sizes FILE` | *(BAM header)* | Tab-delimited chromosome sizes file (`chrom\tsize`). When omitted, chromosome names and sizes are read from the BAM header. |
| `--regions FILE` | *(whole genome)* | BED file of genomic regions to restrict processing to. Only reads overlapping these intervals are examined, which can substantially reduce runtime for targeted analyses. Overlapping intervals within the same chromosome are merged automatically to prevent double-counting. |

### Signal

| Argument | Default | Description |
|---|---|---|
| `--event {edit,tn5}` | `edit` | What the track counts. `edit`: deamination events (C→T or G→A reference mismatches), per fragment. `tn5`: Tn5 insertion sites — the 5′ end of every passing read, shifted +4 bp (forward) or −5 bp from the exclusive end (reverse) onto the centre of the 9-bp duplication. A paired-end fragment therefore contributes **both** of its ends and a single-end read only its start (for a reverse read, the right end of its alignment). `tn5` works with `--mode count` only; `--min_baseq` does not apply to it. See the *Tn5 cut sites* section of [Algorithm](../algorithm.md). |
| `--forward_shift INT` | `+4` | `--event tn5` only. Added to a forward read's alignment start to give its cut site. |
| `--reverse_shift INT` | `−5` | `--event tn5` only. Added to a reverse read's **exclusive** alignment end to give its cut site. With the defaults both reads of one insertion land on the centre of the 9-bp duplication. `--forward_shift 0 --reverse_shift -1` gives raw, unshifted 5′ ends (a reverse read's 5′ base is `end − 1`); `0` and `0` suits a BAM already shifted +4/−5 upstream. Either one set with `--event edit` logs a warning and is ignored. |
| `--extend_size INT` | `0` | Symmetrically extend each detected deamination site by INT base pairs in both directions before writing to the BigWig. A value of 50 means each event at position *p* contributes signal to [*p*−50, *p*+50]. Implemented as a box-kernel convolution, so the signal at a position equals the number of events within `extend_size` bases. |
| `--normalize` | *(off)* | In count mode (the default), scale every value by `scale_factor / (genome-wide total count)` so the written track sums to `--scale_factor` — reads/counts-per-million-style normalization that makes samples comparable regardless of editing depth. With `--extend_size 0` this is counts per `scale_factor` edits. Ignored in `--mode ratio`. |
| `--scale_factor FLOAT` | `1000000` | Target total for `--normalize` (1e6 gives per-million values). |

### Quality filters

| Argument | Default | Description |
|---|---|---|
| `--min_mapq INT` | `20` | Minimum read mapping quality (MAPQ). Reads strictly below this value are skipped entirely. |
| `--min_baseq INT` | `20` | Minimum base quality (phred score) at a candidate position. Individual bases below this value are not counted even if the read passes MAPQ filtering. |

Regardless of these thresholds, the following reads are always excluded:

- Unmapped reads (flag `0x4`)
- PCR/optical duplicates (flag `0x400`)
- QC-failed reads (flag `0x200`)
- Secondary alignments (flag `0x100`)
- Supplementary alignments (flag `0x800`)

### Performance

| Argument | Default | Description |
|---|---|---|
| `--threads INT` | `1` | Number of worker processes. Each region is processed independently, so the work spreads across cores. |

### Global option (before the subcommand)

| Argument | Default | Description |
|---|---|---|
| `--log_level LEVEL` | `INFO` | Logging verbosity. One of `DEBUG`, `INFO`, `WARNING`, `ERROR`. |

## Input file preparation

Before running `bam2bw`, your BAM and FASTA files must be sorted and indexed.

```bash
# Sort and index BAM
samtools sort -o sample.sorted.bam sample.bam
samtools index sample.sorted.bam

# Index reference FASTA
samtools faidx hg38.fa

# (Optional) Generate chromosome sizes file from FASTA index
cut -f1,2 hg38.fa.fai > hg38.chrom.sizes
```

## Examples

### Whole-genome run

```bash
deamtools bam2bw \
    --bam sample.sorted.bam \
    --fasta hg38.fa \
    --out_dir results \
    --out_name sample
```

Writes `results/sample.bw`. Uses default MAPQ ≥ 20 and base quality ≥ 20. Chromosome sizes are read from the BAM header.

### Restrict to peaks and use stricter quality filters

```bash
deamtools bam2bw \
    --bam sample.sorted.bam \
    --fasta hg38.fa \
    --regions peaks.bed \
    --min_mapq 30 \
    --min_baseq 30 \
    --out_dir results \
    --out_name sample_peaks
```

Only reads overlapping intervals in `peaks.bed` are processed, which is much faster than a whole-genome run when peaks cover a small fraction of the genome.

### Extend each deamination event by 50 bp

```bash
deamtools bam2bw \
    --bam sample.sorted.bam \
    --fasta hg38.fa \
    --extend_size 50 \
    --out_dir results \
    --out_name sample_extended
```

Each C→T event contributes signal to a 101-bp window centred on the event position. Useful when the raw per-base signal is too sparse for downstream peak calling.

### Parallel chromosome processing with explicit chromosome sizes

```bash
deamtools bam2bw \
    --bam sample.sorted.bam \
    --fasta hg38.fa \
    --chrom_sizes hg38.chrom.sizes \
    --threads 8 \
    --out_dir results \
    --out_name sample
```

### Enable debug logging

```bash
deamtools --log_level DEBUG bam2bw \
    --bam sample.sorted.bam \
    --fasta hg38.fa \
    --out_dir results \
    --out_name sample
```

Note that `--log_level` is a global flag and must appear **before** the subcommand name.

## Output

The output is a BigWig file in variable-step format. Only positions with non-zero signal are written, keeping file sizes small. The BigWig is directly compatible with genome browsers (IGV, UCSC) and downstream tools (deepTools, MACS3, etc.).

The value at each position is the number of deamination events observed there. If `--extend_size` is used, the value at position *p* is the number of raw events within `extend_size` bases of *p*.

## Choosing parameters

**`--min_mapq`** — A value of 20 (default) retains reads with ≥ 99% mapping accuracy. For repetitive regions or multi-mapping reads, raise to 30. Setting to 0 disables MAPQ filtering.

**`--min_baseq`** — A value of 20 (default) corresponds to 99% base-call accuracy. Lowering increases sensitivity but also increases noise from sequencing errors. Raising above 30 can be overly strict for older sequencing data.

**`--extend_size`** — Set to 0 (default) for the raw single-base signal. For footprinting or broad accessibility analysis, values of 50–200 bp are typical. The optimal value depends on the expected size of accessible regions in your assay.

**`--threads`** — The number of worker processes. Regions are grouped into batches of neighbouring intervals — about eight per worker, each opening the BAM once — so a peak file with tens of thousands of intervals costs a few hundred file opens rather than one per interval. On a whole-genome run each chromosome is its own batch, so each worker holds one chromosome's signal array at a time and memory grows with the worker count. For a human genome, 8–16 is a reasonable range.