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 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.
--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.

# 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

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

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

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

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

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.