Metadata-Version: 2.4
Name: bamnado
Version: 0.9.0
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Rust
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: MacOS :: MacOS X
Classifier: Operating System :: Microsoft :: Windows
Requires-Dist: numpy>=2
Requires-Dist: pytest>=7.0 ; extra == 'test'
Requires-Dist: pytest-cov>=4.0 ; extra == 'test'
Provides-Extra: test
Summary: Python bindings for BamNado
Keywords: bioinformatics,scATAC-seq,MCC,BAM
Home-Page: https://github.com/alsmith151/BamNado
Author-email: Alastair Smith <alastair.smith@ndcls.ox.ac.uk>
License: MIT OR Apache-2.0
Requires-Python: >=3.10
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Documentation, https://alsmith151.github.io/BamNado/
Project-URL: Homepage, https://github.com/alsmith151/BamNado
Project-URL: Repository, https://github.com/alsmith151/BamNado

# BamNado

High-performance tools and utilities for working with BAM and BigWig files in modern genomics workflows. BamNado is written in Rust for speed and low memory use and provides both a command-line interface and Python bindings.

## Overview

BamNado is designed for efficient, streaming manipulation of BAM files and signal tracks. It focuses on fast coverage generation, flexible filtering, and lightweight post-processing of bedGraph and BigWig data.

**Common use cases:**

- Rapid generation of coverage tracks from large BAM files
- Filtering reads by tags or barcodes to produce targeted BigWigs
- Fragment-aware coverage for ATAC-seq and related assays
- BigWig comparison and aggregation across samples
- Post-processing of binned signal tracks for visualization

Useful in workflows including single-cell and Micro-Capture-C (MCC), and many others.

**Documentation:** [API docs on docs.rs](https://docs.rs/bamnado/latest/bamnado/)

## Features

- High-performance, streaming implementations in Rust
- Cross-platform support (Linux, macOS, Windows)
- BAM → bedGraph / BigWig coverage generation
- Fragment-aware and strand-specific pileups
- Read filtering by mapping quality, length, strand, fragment size, tags, and barcodes
- BigWig comparison (subtraction, ratio, log-ratio)
- BigWig aggregation (sum, mean, median, min, max)
- bedGraph post-processing with `collapse-bedgraph`
- Python bindings for selected functionality

## Installation

### Pre-built binaries (recommended)

Download the appropriate binary from the [releases page](https://github.com/alsmith151/BamNado/releases).

After downloading, move to your user bin directory (no root required):

```bash
chmod +x bamnado
mv bamnado ~/.local/bin/bamnado
```

Ensure `~/.local/bin` is on your `PATH` (add to `~/.bashrc` or `~/.zshrc` if needed):

```bash
export PATH="$HOME/.local/bin:$PATH"
```

Alternatively, you can run the binary directly from the download location:

```bash
./bamnado --help
```


### cargo binstall

If you have [`cargo-binstall`](https://github.com/cargo-bins/cargo-binstall) installed, this downloads the pre-built binary automatically:

```bash
cargo binstall bamnado
```

### Docker

```bash
docker pull ghcr.io/alsmith151/bamnado:latest
docker run --rm ghcr.io/alsmith151/bamnado:latest --help
```

Images are available for `linux/amd64` and `linux/arm64`.

### Cargo

If you have Rust installed:

```bash
cargo install bamnado
```

### Build from source

```bash
git clone https://github.com/alsmith151/BamNado.git
cd BamNado
cargo build --release
```

### Optional dependency: samtools

`samtools` is not required but is **strongly recommended** if your BAM files have non-standard or incomplete headers (e.g. files produced by CellRanger). BamNado automatically falls back to `samtools view -H` to parse the header when the built-in parser fails. Without `samtools` on your `PATH`, BamNado will error on such files.

Install via conda or your system package manager:

```bash
conda install -c bioconda samtools
# or
brew install samtools
```

## Command-line interface

### Getting help

```bash
bamnado --help                      # List all commands
bamnado <command> --help            # Help for a specific command
```

### Available commands

| Command | Aliases | Description |
| ------- | ------- | ----------- |
| `bam-coverage` | `coverage` | BAM → bedGraph / BigWig coverage |
| `multi-bam-coverage` | `multi-coverage` | Merge coverage from multiple BAM files into one track |
| `split` | — | Split BAM by tag or barcode |
| `split-exogenous` | `split-spikein` | Split BAM into endogenous / exogenous reads |
| `modify` | — | Filter and/or adjust reads in a BAM file |
| `bigwig-compare` | `compare-bigwigs` | Compare two BigWig files bin by bin |
| `bigwig-aggregate` | `aggregate-bigwigs` | Aggregate multiple BigWig files into one track |
| `bigwig-infer-scale` | `infer-scale` | Infer scaling factor and library size from a normalised BigWig |
| `bam-normalize` | `norm-factors` | Compute between-sample scaling factors from BAM files |
| `collapse-bedgraph` | `collapse` | Collapse adjacent equal-score bins in a bedGraph |

### Read filtering

All coverage commands share common read filter flags:

| Flag | Default | Description |
| ---- | ------- | ----------- |
| `--strand` | `both` | `forward`, `reverse`, or `both` |
| `--proper-pairs` | off | Keep only properly-paired reads |
| `--min-mapq` | 20 | Minimum mapping quality |
| `--min-length` | 20 | Minimum read length (bp) |
| `--max-length` | 1000 | Maximum read length (bp) |
| `--min-fragment-len` | — | Minimum insert size (bp); paired-end only |
| `--max-fragment-len` | — | Maximum insert size (bp); paired-end only |
| `--blacklist` | — | BED file of regions to exclude |
| `--barcode-allowlist` | — | Text file of cell barcodes (one per line) |
| `--read-group` | — | Keep only this read group |
| `--tag` / `--tag-value` | — | Keep reads where TAG == VALUE |
| `--ignore-duplicates` | off | Exclude PCR/optical duplicate reads |
| `--ignore-secondary` | off | Exclude secondary alignments |
| `--ignore-supplementary` | off | Exclude supplementary alignments |

### Coverage-specific flags

| Flag | Default | Description |
| ---- | ------- | ----------- |
| `--normalize` | `raw` | Signal normalization method: `raw`, `rpkm`, or `cpm` |
| `--fragment-counts` | off | Count fragments instead of individual read alignments |
| `--ignore-scaffolds` | off | Skip scaffold or unplaced chromosomes |
| `--threads` | `6` | Threads used when writing BigWig output |

### Examples

#### Generate basic coverage

```bash
bamnado bam-coverage \
  --bam sample.bam \
  --output coverage.bedgraph \
  --bin-size 100
```

#### Generate high-quality, normalized coverage

```bash
bamnado bam-coverage \
  --bam sample.bam \
  --output coverage_hq.bw \
  --bin-size 100 \
  --normalize rpkm \
  --min-mapq 30 \
  --proper-pairs
```

#### Extract nucleosome-free regions from ATAC-seq

```bash
bamnado bam-coverage \
  --bam atac.bam \
  --output nfr_forward.bw \
  --bin-size 10 \
  --strand forward \
  --min-fragment-len 100 \
  --max-fragment-len 200 \
  --fragment-counts \
  --min-mapq 20
```

#### Filter by cell barcode (single-cell)

```bash
bamnado bam-coverage \
  --bam possorted_genome_bam.bam \
  --output cell_coverage.bw \
  --bin-size 100 \
  --barcode-allowlist barcodes.txt \
  --fragment-counts
```

#### Filter by SAM tag (MCC viewpoint)

```bash
bamnado bam-coverage \
  --bam mcc.bam \
  --output BCL2_viewpoint.bw \
  --bin-size 50 \
  --tag VP \
  --tag-value BCL2 \
  --fragment-counts \
  --min-mapq 30
```

#### Compare two BigWig files

```bash
bamnado bigwig-compare \
  --bw1 sample_treated.bw \
  --bw2 sample_control.bw \
  --comparison log-ratio \
  --pseudocount 1e-3 \
  --threads 6 \
  -o treated_vs_control.bw
```

#### Aggregate multiple BigWig files

```bash
bamnado bigwig-aggregate \
  --bigwigs sample1.bw sample2.bw sample3.bw \
  --method mean \
  --threads 6 \
  -o mean_coverage.bw
```

#### Simplify bedGraph file

```bash
bamnado collapse-bedgraph \
  --input signal.bedgraph \
  --output signal.collapsed.bedgraph
```

#### Normalize between ChIP-seq samples

```bash
bamnado bam-normalize \
  --bams chip1.bam chip2.bam chip3.bam \
  --method csaw-background \
  --bin-size 10000 \
  --blacklist blacklist.bed
```

Apply the reported `scale_factor` for each sample to `bam-coverage -f` (or
`bigwig-aggregate --scale-factors`) when generating tracks for that sample.

## Between-sample normalization

Comparing ChIP-seq / CUT&Tag / CUT&RUN samples side by side requires removing
composition and technical bias, not just sequencing depth. A sample with a stronger
pulldown has more of its library consumed by peaks, which deflates its background and
makes every other region look artificially low relative to a shallower or weaker-pulldown
sample. Simple total-count (CPM) scaling does not correct for this — it's the same
problem `csaw::normFactors` / `edgeR::calcNormFactors` (TMM) solve in R/Bioconductor.

`bamnado bam-normalize` bins the genome, counts filtered reads (or fragments) per bin per
sample, and reports a per-sample scaling factor plus the underlying bin count matrix for
downstream sanity checks (e.g. M-A plots). It reports factors only — it does not generate
bigwigs itself.

### Method choice

| Method | When to use |
| ------ | ----------- |
| `csaw-background` (default) | Typical composition bias: TMM restricted to large background windows with the most-enriched bins excluded |
| `tmm` | TMM over all bins, without excluding enriched regions |
| `cpm` | Naive depth-only baseline, useful for comparison against the other methods |
| `median-of-ratios` | A robust DESeq2-style alternative estimator |
| `spike-in` | When a spike-in genome was used and is trusted more than a background estimate |

### Applying the result

`scale_factor` is centred on 1 (the geometric mean across samples), so it can be passed
directly to `bam-coverage -f <factor>` or `bigwig-aggregate --scale-factors` for each
corresponding sample.

Use `--counts-out counts.tsv` to write the underlying bin count matrix (chrom, start, end,
one column per sample) for QC, e.g. plotting an M-A plot to confirm the fitted factor sits
at the centre of the trimmed M-value cloud.

### Spike-in caveat

The `spike-in` method reads mapped-read counts directly from each BAM's index (no BAM
record scan), so those counts are **unfiltered** by MAPQ, duplicates, or
secondary/supplementary status — the same basis as `bam-coverage`'s internal mapped-read
count. `--exogenous-prefix` is required and `--bin-size` is ignored for this method.

## Python API

Installing the Python package also installs the `bamnado` CLI, including all Rust CLI
commands and flags:

```bash
pip install bamnado
bamnado bam-coverage --help
```

BamNado provides Python bindings for high-performance BAM signal generation with flexible read filtering. Install from PyPI:

```bash
pip install bamnado
# or
uv pip install bamnado
```

### Quick start

```python
import bamnado
import numpy as np

# Generate coverage from a BAM file
signal = bamnado.get_signal_for_chromosome(
    bam_path="input.bam",
    chromosome_name="chr1",
    bin_size=50,
)
print(f"Mean coverage: {np.mean(signal):.2f}")
```

### ReadFilter class

Customize read filtering with the `ReadFilter` class. All parameters are optional:

```python
import bamnado

# Create a filter for high-quality, properly-paired reads
strict_filter = bamnado.ReadFilter(
    min_mapq=30,
    proper_pair=True,
    min_length=50,
)

# Apply the filter
signal = bamnado.get_signal_for_chromosome(
    bam_path="input.bam",
    chromosome_name="chr1",
    bin_size=100,
    read_filter=strict_filter,
)
```

#### Filter parameters

| Parameter | Type | Default | Description |
| --------- | ---- | ------- | ----------- |
| `min_mapq` | `int` | `0` | Minimum mapping quality |
| `proper_pair` | `bool` | `False` | Require properly paired reads |
| `min_length` | `int` | `0` | Minimum read length (bp) |
| `max_length` | `int` | `1000` | Maximum read length (bp) |
| `strand` | `str` | `"both"` | `"forward"`, `"reverse"`, or `"both"` |
| `min_fragment_length` | `int \| None` | `None` | Minimum insert size (bp); paired-end only |
| `max_fragment_length` | `int \| None` | `None` | Maximum insert size (bp); paired-end only |
| `blacklist_bed` | `str \| None` | `None` | BED file of excluded regions |
| `whitelisted_barcodes` | `list[str] \| None` | `None` | Cell barcodes to include (CB tag) |
| `read_group` | `str \| None` | `None` | Read group to keep (RG tag) |
| `filter_tag` | `str \| None` | `None` | SAM tag to filter on (e.g. `"VP"`) |
| `filter_tag_value` | `str \| None` | `None` | Required value for `filter_tag` |

### Usage examples

#### Basic coverage

```python
import bamnado
import numpy as np

signal = bamnado.get_signal_for_chromosome(
    bam_path="sample.bam",
    chromosome_name="chr1",
    bin_size=50,
    scale_factor=1.0,
    use_fragment=False,
)
print(f"Signal shape: {signal.shape}, type: {signal.dtype}")
print(f"Mean: {np.mean(signal):.2f}, Max: {np.max(signal):.2f}")
```

#### Nucleosome-free regions (ATAC-seq)

```python
import bamnado

# Forward-strand reads, 100–200 bp fragments
nfr_filter = bamnado.ReadFilter(
    strand="forward",
    min_fragment_length=100,
    max_fragment_length=200,
    min_mapq=20,
)

nfr_signal = bamnado.get_signal_for_chromosome(
    bam_path="atac.bam",
    chromosome_name="chr1",
    bin_size=10,
    use_fragment=True,
    read_filter=nfr_filter,
)
```

#### Barcode-filtered coverage (single-cell)

```python
import bamnado

# Get coverage for specific cell barcodes
cell_barcodes = ["ACGTACGT-1", "TGCATGCA-1", "AAATAAAA-1"]

barcode_filter = bamnado.ReadFilter(
    whitelisted_barcodes=cell_barcodes,
    min_mapq=25,
)

signal = bamnado.get_signal_for_chromosome(
    bam_path="possorted_genome_bam.bam",
    chromosome_name="chr1",
    bin_size=100,
    use_fragment=True,
    read_filter=barcode_filter,
)
```

#### Tag-filtered coverage (e.g. MCC viewpoint)

```python
import bamnado

# Get coverage for specific viewpoint
vp_filter = bamnado.ReadFilter(
    filter_tag="VP",
    filter_tag_value="BCL2",
    min_mapq=30,
)

vp_signal = bamnado.get_signal_for_chromosome(
    bam_path="mcc.bam",
    chromosome_name="chr1",
    bin_size=50,
    use_fragment=True,
    read_filter=vp_filter,
)
```

#### High-confidence coverage

```python
import bamnado

# High-quality, properly-paired reads
confident_filter = bamnado.ReadFilter(
    min_mapq=30,
    proper_pair=True,
    min_length=50,
    min_fragment_length=100,
    max_fragment_length=500,
)

signal = bamnado.get_signal_for_chromosome(
    bam_path="sample.bam",
    chromosome_name="chr1",
    bin_size=50,
    scale_factor=1e6,  # CPM normalization
    use_fragment=True,
    read_filter=confident_filter,
)
```

### get_signal_for_chromosome parameters

| Parameter | Type | Default | Description |
| --------- | ---- | ------- | ----------- |
| `bam_path` | `str` | — | Path to indexed BAM file |
| `chromosome_name` | `str` | — | Chromosome to process (e.g. `"chr1"`) |
| `bin_size` | `int` | `50` | Bin width in base pairs |
| `scale_factor` | `float` | `1.0` | Linear scaling factor |
| `use_fragment` | `bool` | `False` | Use fragment (pair) coordinates instead of read |
| `ignore_scaffold_chromosomes` | `bool` | `True` | Skip non-standard chromosomes |
| `read_filter` | `ReadFilter \| None` | `None` | Optional read filter |

**Returns:** `numpy.ndarray` of dtype `float32` with length = chromosome size / bin_size

**Note:** Fragment length filtering (`min_fragment_length`, `max_fragment_length`) requires paired-end BAM files and will raise `ValueError` on single-end data.

### Between-sample scaling factors

`compute_scale_factors` mirrors the `bam-normalize` CLI command — see
[Between-sample normalization](#between-sample-normalization) for the underlying method.

```python
import bamnado

result = bamnado.compute_scale_factors(
    ["chip1.bam", "chip2.bam", "chip3.bam"],
    method="csaw-background",
    bin_size=10000,
)
print(result.sample_names, result.scale_factors)

# Include the bin count matrix for QC (e.g. an M-A plot)
result = bamnado.compute_scale_factors(
    ["chip1.bam", "chip2.bam"],
    method="tmm",
    return_counts=True,
)
counts = result.counts()  # (n_bins, n_samples) uint64 array, or None
```

#### compute_scale_factors parameters

| Parameter | Type | Default | Description |
| --------- | ---- | ------- | ----------- |
| `bam_paths` | `list[str]` | — | Input BAM file paths, one per sample |
| `method` | `str` | `"csaw-background"` | `"tmm"`, `"csaw-background"`, `"cpm"`, `"median-of-ratios"`, or `"spike-in"` |
| `bin_size` | `int` | `10000` | Genomic bin width in bp; ignored for `"spike-in"` |
| `sample_names` | `list[str] \| None` | `None` | Defaults to each BAM's file stem, disambiguated if duplicated |
| `exclude_top_percent` | `float` | `5.0` | `"csaw-background"` only: percent of highest-mean bins dropped before estimating |
| `reference_sample` | `str \| None` | `None` | TMM reference sample; defaults to the one closest to the geometric mean |
| `logratio_trim` | `float` | `0.3` | TMM trim fraction for M-values |
| `sum_trim` | `float` | `0.05` | TMM trim fraction for A-values |
| `exogenous_prefix` | `str \| None` | `None` | Required for `method="spike-in"` |
| `use_fragment` | `bool` | `False` | Count fragments (pairs) instead of individual reads |
| `ignore_scaffold_chromosomes` | `bool` | `True` | Skip scaffold / unplaced chromosomes |
| `read_filter` | `ReadFilter \| None` | `None` | Applied per sample |
| `return_counts` | `bool` | `False` | Populate `ScaleFactorResult.counts()` |

**Returns:** `ScaleFactorResult` with `method`, `sample_names`, `library_sizes`,
`norm_factors`, `scale_factors`, `reference_sample`, `n_bins_total`, `n_bins_used`, and
(when `return_counts=True`) `counts()`, `bin_chroms`, `bin_starts`, `bin_ends`.

## Development

```bash
cargo build --release
cargo test
```

## License

Apache-2.0 OR MIT

