Metadata-Version: 2.4
Name: alignable
Version: 0.1.0
Summary: Aligner-specific mappability and splice junction artifact tracks for RNA-seq reference genomes
Keywords: bioinformatics,genomics,rna-seq,mappability,splice-junction,alignment
Author-Email: Matthew Iyer <mkiyer@umich.edu>
Maintainer-Email: Matthew Iyer <mkiyer@umich.edu>
License-Expression: MIT
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: MacOS
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: C++
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: Implementation :: CPython
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Project-URL: Homepage, https://github.com/mkiyer/alignable
Project-URL: Repository, https://github.com/mkiyer/alignable.git
Project-URL: Issues, https://github.com/mkiyer/alignable/issues
Project-URL: Changelog, https://github.com/mkiyer/alignable/blob/main/CHANGELOG.md
Requires-Python: >=3.11
Requires-Dist: pysam>=0.22
Requires-Dist: zarr>=3.0
Requires-Dist: numpy>=1.24
Requires-Dist: pyarrow>=14
Requires-Dist: click>=8
Requires-Dist: pyyaml>=6
Provides-Extra: viz
Requires-Dist: pybigtools>=0.2.5; extra == "viz"
Requires-Dist: matplotlib; extra == "viz"
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Description-Content-Type: text/markdown

# alignable

Aligner-specific mappability and splice junction artifact data for
reference genomes. Generates per-base mappability tracks and splice
junction blacklists using synthetic paired-end reads aligned with
STAR, minimap2, or hisat2.

## Overview

alignable tiles a reference genome with synthetic paired-end fragments,
aligns them back to the genome, and scores each position for
mappability. The output is a Zarr v3 store with four channels per
position per read-length bin:

| Channel | Description |
|---------|-------------|
| `n_frag` | Total fragments covering this position |
| `n_frag_ontarget` | Fragments with both mates mapping to correct position |
| `nh_ontarget` | Sum of NH (number of hits) for on-target fragments |
| `n_spliced` | Fragments with splice junction (CIGAR N-op) artifacts |

From these channels, two key metrics are derived:

- **Binary mappability** = `n_frag_ontarget / n_frag` — fraction of
  fragments that map uniquely to the correct position
- **Fractional mappability** = `n_frag_ontarget / nh_ontarget` — inverse
  of the effective copy number (≈ 1/CN for duplicated regions)

## Breaking change notice — store format

This release of `alignable` **does not read stores produced by older
versions**. Three things changed and there is no automatic
back-compat path:

1. **Splice blacklist on disk.** The aggregated blacklist is now a
   single zstd-compressed Feather v2 file at
   `<store>/mappability.zarr/splice_blacklist.feather`, replacing the
   former `splice_blacklist/` Zarr group (per-column arrays + JSON
   strings in `.zattrs`). New layout is ~5× smaller and ~17× faster
   to load. The legacy Zarr-group reader has been removed.
2. **Per-chunk splice artifacts.** Only `chunk_*.feather` files are
   read from `splice/`. The previous `chunk_*.parquet` and
   `chunk_*.json[.gz]` formats are no longer recognized.
3. **Packaged zip location.** `alignable package` now writes the
   `.zarr.zip` as a **sibling** of the output directory
   (`<parent>/<output_dir_name>.zarr.zip`), not inside it. The
   directory store and the zip are independent — there is no
   silent fallback between them when you call `alignable.open()`.

**To migrate an existing run:** if you still have the per-chunk
`splice/chunk_*.feather` files, re-run `alignable finalize` (or
`scripts/reaggregate_splice.py`) to rebuild the blacklist in the
new format. If you only have an old `splice_blacklist/` Zarr group
and the `splice/` chunks are gone, you must rebuild the store with
`alignable compute`.

## Installation

Alignable is distributed as source. You install it into a conda
environment that supplies the C++ toolchain, htslib, Python
dependencies, and (optionally) the aligner.

### Quickstart (read-only / query use)

If you only need to **query an existing store** produced by someone
else — i.e. you are not running the alignment pipeline — the
install is minimal:

```bash
# Clone the repository
git clone https://github.com/<your-org>/alignable.git
cd alignable

# Create the env (provides Python 3.13, numpy, zarr 3, pysam, etc.)
mamba env create -f mamba_env.yaml
conda activate alignable

# Install alignable (builds the C++ extensions; required for the package
# import to work, even if you only call the read-only API)
pip install --no-build-isolation -ve .

# Verify
python -c "import alignable; print(alignable.__version__)"
```

No aligner is needed for query-only use.

### Full install (compute pipeline)

For running the pipeline (`alignable compute` etc.) you also need an
aligner:

```bash
# After the steps above:
mamba install -n alignable -c bioconda star      # recommended for RNA-seq
# or
mamba install -n alignable -c bioconda hisat2
# minimap2 is already included in mamba_env.yaml
```

### What `pip install` builds

Two nanobind C++ extensions are compiled in place using
scikit-build-core:

- `_tiler_fast` — batch FASTQ generation (xoshiro256** PRNG, ~33× speedup)
- `_scorer_fast` — htslib-based BAM scoring (~250× speedup)

Both extensions build automatically using the compiler (`cxx-compiler`)
and `htslib` from the conda environment. You do **not** need
system-level htslib, GCC, or `module load`.

### Requirements summary

| Component | Source | Notes |
|-----------|--------|-------|
| Python ≥ 3.10 (3.13 in env) | conda | |
| C++ compiler (GCC ≥ 12) | conda (`cxx-compiler`) | For nanobind extensions |
| cmake ≥ 3.15 | conda | Build system |
| htslib ≥ 1.19 | conda | For `_scorer_fast` C++ extension |
| zarr ≥ 3.0 | conda | Store format |
| pysam, numpy, click, pyyaml | conda | Core deps |
| STAR / minimap2 / hisat2 | conda (bioconda) | Pipeline only — not needed for queries |
| pybigtools | pip | Optional, only for `export-bigwig` |

### Requirements summary

| Component | Source | Notes |
|-----------|--------|-------|
| Python ≥ 3.10 | conda | |
| C++ compiler (GCC ≥ 12) | conda (`cxx-compiler`) | For nanobind extensions |
| cmake ≥ 3.15 | conda | Build system |
| htslib ≥ 1.19 | conda | For `_scorer_fast` C++ extension |
| pysam ≥ 0.22 | conda | BAM I/O |
| zarr ≥ 3.0 | conda | Output store format |
| numpy, pyarrow, click | conda | Core dependencies |
| pyyaml ≥ 6 | conda | Config file support |
| STAR / minimap2 / hisat2 | conda (bioconda) | At least one required |

## Quick start

### Single-node genome-wide run (STAR)

```bash
alignable compute \
  --genome /path/to/genome.fa \
  --aligner star \
  --aligner-index /path/to/star_index \
  --aligner-config star_config.txt \
  --read-lengths 50,75,100,125,150,200 \
  --threads 24 \
  -j 6 \
  --max-chunk-size 10000000 \
  --output-dir /path/to/output
```

When `-j > 1` with STAR, alignable automatically loads the genome index
into shared memory (`--genomeLoad LoadAndKeep`) so all parallel STAR
instances share a single copy (~30 GB), and removes it when finished.

### Using a config file

Instead of passing all parameters on the command line, use a YAML
config file. The config is saved to `config.yaml` in the output
directory, serving as a record of exactly how the run was produced.

```yaml
# alignable_config.yaml
genome: /path/to/genome.fa
aligner: star
aligner-index: /path/to/star_index
aligner-config: star_config.txt
read-lengths: "50,75,100,150"
frag-len-mean: 300
frag-len-sd: 50
frag-len-min: 1
frag-len-max: 1000
error-rate: 0.01
coverage: 1
tolerance: 5
threads: 24
parallel: 6
max-chunk-size: 10000000
output-dir: /path/to/output
seed: 42
```

```bash
alignable --config alignable_config.yaml compute
```

CLI flags override config file values:

```bash
# Use config but override threads and parallelism
alignable --config alignable_config.yaml compute --threads 48 -j 12
```

## Parameters

### Fragment generation

alignable generates synthetic paired-end fragments by tiling the
reference genome. Each fragment consists of two reads (R1 forward,
R2 reverse-complement) extracted from a genomic interval. Fragment
sizes are drawn from a truncated Gaussian distribution.

| Parameter | CLI flag | Default | Description |
|-----------|----------|---------|-------------|
| Read lengths | `--read-lengths` | *(required)* | Comma-separated list of read lengths to simulate (e.g. `75,150`). Each length is a separate "bin" — the tool generates a complete set of fragments for every bin. |
| Fragment length mean | `--frag-len-mean` | 300 | Mean fragment length in bp. This is the total genomic span of the fragment, not the sequenced portion. Typical values: 200–400 for standard RNA-seq libraries. |
| Fragment length SD | `--frag-len-sd` | 50 | Standard deviation of fragment lengths in bp. Fragment sizes are drawn from a Gaussian clamped to `[frag_len_min, min(frag_len_max, frag_len_mean + 3 × frag_len_sd)]`. |
| Fragment length min | `--frag-len-min` | 1 | Hard floor on fragment length. Fragments shorter than this are clamped up. |
| Fragment length max | `--frag-len-max` | 1000 | Hard ceiling on fragment length. Prevents excessively long fragments when SD is large. |
| Coverage | `--coverage` | 1 | Number of independent fragments generated per genomic position per read-length bin. Higher values reduce stochastic noise in mappability estimates but increase runtime proportionally. Coverage of 1 is sufficient for most applications since mappability is primarily determined by sequence uniqueness. |
| Error rate | `--error-rate` | 0.01 | Per-base substitution error rate applied to synthetic reads. Simulates sequencing errors. Set to 0 for error-free reads. |

### Scoring

| Parameter | CLI flag | Default | Description |
|-----------|----------|---------|-------------|
| Tolerance | `--tolerance` | 5 | Maximum allowed distance (in bp) between the expected and observed alignment position for a fragment to be considered "on-target". Accounts for soft-clipping and small indels in the alignment. |

### Aligner

| Parameter | CLI flag | Default | Description |
|-----------|----------|---------|-------------|
| Aligner | `--aligner` | *(required)* | Alignment program: `star`, `minimap2`, or `hisat2`. |
| Aligner index | `--aligner-index` | *(required)* | Path to the pre-built aligner index. |
| Aligner config | `--aligner-config` | None | Plain-text file with aligner-specific arguments (one per line). See [Aligner config files](#aligner-config-files) below. |

### Execution

| Parameter | CLI flag | Default | Description |
|-----------|----------|---------|-------------|
| Threads | `--threads` | 1 | Total thread budget. When using `-j N`, threads are divided evenly: each worker gets `threads // N` aligner threads. |
| Parallel | `-j`, `--parallel` | 1 | Number of chunks to process concurrently. Each worker runs its own aligner subprocess. For STAR with `-j > 1`, shared memory is used automatically. |
| Max chunk size | `--max-chunk-size` | None | Split chromosomes longer than this into multiple chunks (bp). Enables parallelism across large chromosomes. 10 Mbp is a good default. |
| Seed | `--seed` | None | RNG seed for reproducible fragment generation. |
| Output dir | `--output-dir` | *(required)* | Output directory. The Zarr store, chunk plan, and config are written here. |

### Fragment size distribution explained

The fragment size distribution controls how genomic intervals
are sampled:

```
                    frag_len_mean = 300
                        │
     ┌──────────────────┼──────────────────┐
     │                  │                  │
     ▼                  ▼                  ▼
  ───┤──────────────────────────────────┤───
   frag_len_mean          │          frag_len_mean
   - 3 × frag_len_sd      │          + 3 × frag_len_sd
   (= 150)              │          (= 450)
                        │
          Fragment: [────R1────]........[────R2────]
                    read_length         read_length

  Minimum fragment length = frag_len_min (default: 1)
  Maximum fragment length = min(frag_len_max, frag_len_mean + 3 × frag_len_sd)
```

For a standard RNA-seq library with 150 bp reads:
- `--frag-len-mean 300 --frag-len-sd 50` → fragments range up to 450 bp
- `--frag-len-mean 200 --frag-len-sd 30` → fragments range up to 290 bp

## Aligner config files

Create a plain-text file with one aligner argument per line. alignable
controls I/O, threading, and output format — only set alignment
parameters:

### STAR config example

```
# star_config.txt — STAR parameters for alignable
--outFilterMultimapNmax 100
--outSAMmultNmax -1
--winAnchorMultimapNmax 100
--chimSegmentMin 10
--chimOutType WithinBAM HardClip
--chimJunctionOverhangMin 10
--chimScoreDropMax 30
--chimScoreJunctionNonGTAG 0
--chimScoreSeparation 1
--chimSegmentReadGapMax 3
--chimMultimapNmax 100
```

**Forbidden STAR parameters** (controlled by alignable):
`--readFilesIn`, `--outSAMtype`, `--outStd`, `--outSAMunmapped`,
`--runThreadN`, `--genomeDir`, `--genomeLoad`, `--outTmpDir`,
`--parametersFiles`

### minimap2 config example

Use the `splice:sr` preset with a junctions BED file for RNA-seq:

```
# mm2_config.txt
-x splice:sr
-j /path/to/junctions.bed
--secondary=yes
```

Generate the junctions BED from a GTF annotation:

```bash
paftools.js gff2bed genes.gtf.gz > junctions.bed
minimap2 -x splice:sr -d genome_splice_sr.mmi genome.fa
```

## Resource planning

### Thread and parallelism allocation

The `--threads` flag sets the **total** thread budget. With `-j N`, each
parallel chunk gets `threads // N` aligner threads.

| Cores | `-j` | Threads/chunk | Notes |
|-------|------|---------------|-------|
| 24 | 1 | 24 | Single chunk, all threads to aligner |
| 24 | 6 | 4 | 6 parallel STAR, shared memory saves ~150 GB |
| 24 | 12 | 2 | 12 parallel, STAR still effective at 2 threads |
| 48 | 8 | 6 | Large node |

### Memory

| Component | RAM |
|-----------|-----|
| STAR genome index (human) | ~30 GB |
| Per STAR worker (overhead) | ~2–3 GB |
| Python tiler + scorer per worker | ~1–2 GB |
| **Total for 6 parallel (shared memory)** | **~50 GB** |
| **Total for 6 parallel (no shared memory)** | **~200 GB** |

STAR shared memory is **critical** for parallel runs. Without it, each
STAR instance loads its own copy of the ~30 GB genome index. When
`-j > 1` with STAR, shared memory is used automatically.

### Chunk sizing

Use `--max-chunk-size` to split large chromosomes for parallelism.
10 Mbp chunks work well — the human genome produces ~580 chunks,
keeping all workers busy. Without `--max-chunk-size`, each chromosome
is one chunk (286 chunks for the human genome).

### Read-length bins

Each read-length bin multiplies the number of synthetic fragments.

| Use case | `--read-lengths` | Fragments/position |
|----------|------------------|--------------------|
| Single length | `75` | 1× |
| Short + long | `75,150` | 2× |
| Comprehensive | `50,75,100,125,150,200` | 6× |

More bins = more accurate effective length correction but longer runtime.

## HPC / SLURM usage

### Three-step workflow for array jobs

```bash
# 1. Initialize store and chunk plan (fast, run on login node)
alignable init \
  --genome genome.fa \
  --aligner star \
  --aligner-config star_config.txt \
  --read-lengths 50,75,100,125,150,200 \
  --max-chunk-size 10000000 \
  --output-dir /scratch/alignable_output

# 2. Submit array job (one task per chunk)
N_CHUNKS=$(python -c "
from alignable.chunks import load_chunk_plan
print(len(load_chunk_plan('/scratch/alignable_output')))
")

sbatch --array=0-$((N_CHUNKS-1)) <<'EOF'
#!/bin/bash
#SBATCH --cpus-per-task=4
#SBATCH --mem=40G
#SBATCH --time=2:00:00

alignable compute-chunk \
  --output-dir /scratch/alignable_output \
  --chunk-id $SLURM_ARRAY_TASK_ID \
  --aligner-index /path/to/star_index \
  --threads 4 \
  --seed 42
EOF

# 3. Finalize (after all array tasks complete)
alignable finalize --output-dir /scratch/alignable_output
```

### Single-node with shared STAR memory

For a single node with many cores, use `alignable compute` directly.
STAR shared memory is handled automatically when `-j > 1`:

```bash
#!/bin/bash
#SBATCH --cpus-per-task=24
#SBATCH --mem=120G

alignable --config config.yaml compute
```

Or equivalently, with the `--genome-load LoadAndKeep` flag for manual
control in array jobs:

```bash
#!/bin/bash
#SBATCH --cpus-per-task=24
#SBATCH --mem=120G

# Load genome into shared memory (once)
python -c "
from alignable.aligner_config import star_load_genome
star_load_genome('/path/to/star_index')
"

# Run chunks with shared memory
for CHUNK_ID in $(seq 0 $((N_CHUNKS-1))); do
  alignable compute-chunk \
    --output-dir /scratch/output \
    --chunk-id $CHUNK_ID \
    --aligner-index /path/to/star_index \
    --genome-load LoadAndKeep \
    --threads 24 \
    --seed 42
done

# Remove genome from shared memory
python -c "
from alignable.aligner_config import star_remove_genome
star_remove_genome('/path/to/star_index')
"
```

## Querying results

### CLI

```bash
# Summary statistics
alignable summary --store /path/to/output

# Binary mappability for a region (BedGraph to stdout)
alignable query \
  --store /path/to/output \
  --region chr22:16000000-50000000 \
  --read-length 75 \
  --metric binary

# Export BigWig for genome browser
alignable export-bigwig \
  --store /path/to/output \
  --read-length 75 \
  --metric binary \
  -o mappability_rl75.bw
```

### Python API

The high-level read-only API lives in `alignable.api`. You can either
import `AlignableStore` directly or use the `alignable.open()`
shortcut.

#### Opening a store

```python
import alignable

# Unpackaged form: pass either the alignable output directory ...
store = alignable.open("/path/to/output_dir")

# ... or the inner directory store ...
store = alignable.open("/path/to/output_dir/mappability.zarr")

# Packaged form: pass the .zarr.zip path explicitly.
store = alignable.open("/path/to/output_dir.zarr.zip")
```

The unpackaged directory store and the packaged `.zarr.zip` are two
**independent** representations of the same data. `alignable.open()`
opens whichever path you give it; there is no silent fallback
between them. By convention `alignable package` writes the zip as a
sibling of the output directory (`<parent>/<output_dir_name>.zarr.zip`),
not inside it.

#### Inspecting metadata

```python
store.aligner            # 'star' / 'minimap2' / 'hisat2'
store.aligner_version    # e.g. '2.7.11b'
store.read_length_bins   # e.g. [50, 75, 100, 125]
store.chromosomes        # {'chr1': 248956422, 'chr2': 242193529, ...}
store.frag_len_mean      # 300
store.error_rate         # 0.01
store.status             # 'complete'
```

#### Per-base data

All region queries are **0-based, half-open** (BED-style) and require
a `read_length` matching one of `store.read_length_bins`.

```python
# Raw 4-channel uint32 counts: shape (end - start, 4)
# Channels: n_frag, n_frag_ontarget, nh_ontarget, n_spliced
counts = store.counts("chr22", 16_000_000, 17_000_000, read_length=75)

# Per-base binary mappability (float32 in [0, 1])
bm = store.binary_mappability("chr22", 16_000_000, 17_000_000, 75)

# Per-base fractional mappability (float32 in [0, 1])
fm = store.fractional_mappability("chr22", 16_000_000, 17_000_000, 75)

# Sum of binary mappability over a region (single float)
ml = store.mappable_length("chr22", 16_000_000, 17_000_000, 75)
```

#### Mappable / unmappable BED intervals

```python
# Runs of mappable bases
regions = store.mappable_regions(
    "chr22", 0, 50_818_468, read_length=75,
    threshold=0.99,        # min binary mappability
    unique=True,           # also require fractional mappability
    unique_threshold=0.99, # min fractional mappability
)
# -> [('chr22', 16050000, 16050873), ('chr22', 16051002, ...), ...]

# Runs of unmappable bases (positions with no on-target fragments,
# or binary mappability ≤ threshold)
unmap = store.unmappable_regions("chr22", 0, 50_818_468, 75, threshold=0.0)
```

#### Splice junction blacklist

The aggregated splice-artifact blacklist is stored as a single
zstd-compressed Feather v2 file inside the Zarr store at
`<store>/mappability.zarr/splice_blacklist.feather`. It is included
automatically when the store is packaged with `alignable package`.

Two read APIs are exposed:

```python
# Zero-copy Arrow table — preferred for large blacklists (~35 M rows
# for a full human genome). dict-encoded `chrom`/`strand`, int32
# numerics; loads in ~1 s from a directory store.
table = store.splice_blacklist_table()
# pyarrow.Table with columns:
#   chrom, intron_start, intron_end, strand,
#   read_length, count, max_anchor_left, max_anchor_right

# Convenience: same data as a list of dicts (slower; materializes
# every row in Python).
bl = store.splice_blacklist()
# [{'chrom': 'chr1', 'intron_start': 12227, 'intron_end': 12612,
#   'strand': '+', 'read_length': 75, 'count': 4,
#   'max_anchor_left': 38, 'max_anchor_right': 38}, ...]
```

To export the blacklist as a gzipped BED file (not produced by
default), use `alignable export-splice-bed --store <path>`.

### Performance: caching and batch queries

`AlignableStore` decompresses Zarr chunks on read. Each chunk covers
1 000 000 genomic positions × 4 read-length bins × 4 channels and is
~64 MB uncompressed, so a single small query (even 100 bp) costs
**~100 ms** of zstd decompression on first access. To make repeated
queries fast, the store keeps an LRU cache of decompressed chunks:

```python
# Default: cache up to 64 chunks (~4 GB RAM) — fits a full human
# chromosome's worth of chunks, so random access within a single
# chromosome becomes effectively cache-resident after a warmup pass.
store = alignable.open("/path/to/output")

# Cap memory: cache up to 4 chunks (~256 MB)
store = alignable.open("/path/to/output", cache_size=4)

# Disable caching
store = alignable.open("/path/to/output", cache_size=0)

# Inspect / clear
store.cache_info()   # CacheInfo(hits=..., misses=..., maxsize=64, currsize=...)
store.clear_cache()  # free memory
```

**Best practices for thousands of small queries:**

1. **Use `batch_counts()`** — it sorts your regions by genomic
   position so each chunk is decompressed at most once:

   ```python
   regions = [("chr1", s, s + 1000) for s in starts]  # thousands of items
   counts_list = store.batch_counts(regions, read_length=75)
   # counts_list[i] corresponds to regions[i] (original order preserved)
   ```

2. **Reuse a single `AlignableStore` instance** across all queries —
   each instance has its own cache.

3. **Sort BED-driven workflows by `(chrom, start)`** before iterating
   — this turns random access into a near-sequential cache-friendly
   pattern even without `batch_counts`.

4. **Size the cache to your working set.** For genome-scale scans
   processing one chromosome at a time, `cache_size=1` is enough.
   For random sampling across many regions of a single chromosome,
   `cache_size` should be ≥ the number of chunks you touch (~50 for
   chr22, ~250 for chr1, ~3 100 for the full human genome).

**Typical performance** (chr22, ~50 Mbp, 51 chunks):

| Workload | No cache | With cache (size=64) |
|---|---|---|
| Single 1 kb query (cold) | 105 ms | 105 ms |
| Single 1 kb query (warm) | 105 ms | < 0.1 ms |
| 1 000 random 1 kb queries | ~130 s | ~7 s (`batch_counts` or sorted access) |
| Full chr22 `binary_mappability` | ~8 s | ~8 s |

## Output structure

A completed `alignable compute` run produces the following layout:

```
<parent>/
  output_dir/
    config.yaml                 # Effective parameters (YAML; reusable as --config)
    chunk_plan.json             # Genomic chunks used for computation
    aligner_config.txt          # Copy of the user's aligner config
    mappability.zarr/           # Zarr v3 directory store (canonical data)
      zarr.json                 # Group metadata
      chr1/counts/              # uint32 [chrom_len, n_rl_bins, 4]
      chr2/counts/
      ...
      splice_blacklist.feather  # Aggregated splice junction blacklist (Feather v2 + zstd)
    splice/
      chunk_0.feather           # Per-chunk raw splice artifacts (Feather only)
      chunk_1.feather
      ...
    chunks_done/                # Per-chunk completion markers
      chrN__start__end.done

  output_dir.zarr.zip           # Optional: created by `alignable package`,
                                # written as a SIBLING of output_dir/
```

`config.yaml` is a complete record of the parameters used and can be
passed back to `--config` for subsequent runs with the same settings.

### What you actually need to keep / transfer

`mappability.zarr/` (unpackaged) and `output_dir.zarr.zip` (packaged)
contain the **same data** — the zip is a `ZIP_STORED` snapshot of the
directory tree, including the embedded `splice_blacklist.feather`. They
live in **different locations** (the zip is a sibling of the output
directory, not inside it), so you never need to delete one to remove
duplication.

| File / dir | Size (human) | Required for queries? | Purpose |
|---|---|---|---|
| `output_dir/mappability.zarr/` | ~27 GB | **Yes** (or the zip) | Canonical directory store. Use this in place. |
| `output_dir.zarr.zip` | ~27 GB | **Yes** (or the dir) | Single-file packaging. Best for transfer. |
| `output_dir/config.yaml` | small | No | Provenance — records exactly how the store was built. Keep it. |
| `output_dir/chunk_plan.json` | small | No | Only needed to resume / extend a partial run. |
| `output_dir/aligner_config.txt` | small | No | Provenance for the aligner parameters. |
| `output_dir/chunks_done/` | small | No | Per-chunk completion markers; only needed to resume a partial run. |
| `output_dir/splice/` | ~16 GB | No | Raw per-chunk Feather artifacts. Only needed to **re-aggregate** the splice blacklist with different parameters. Safe to delete after `finalize`. |

#### Recommended transfer sets

**To use the store on another machine (read-only / queries)** you
need only **one** of the two stores. Pick whichever is more
convenient:

- **Single-file (recommended for transfer):**
  ```
  output_dir.zarr.zip
  output_dir/config.yaml          # optional but recommended for provenance
  ```
  Open with `alignable.open("output_dir.zarr.zip")`. Slightly slower
  to first-open (~1 s vs instant) but trivial to copy with
  `scp`/`rsync`/`s3 cp`.

- **Directory store (best for in-place workloads):**
  ```
  output_dir/mappability.zarr/
  output_dir/config.yaml          # optional
  ```
  Open with `alignable.open("output_dir/mappability.zarr")` or
  `alignable.open("output_dir")`. Faster random access on local disk
  but contains thousands of small files — use `tar`/`rsync -a` for
  transfer.

**To re-aggregate the splice blacklist later** (e.g. with a different
`--min-count`), additionally keep:

```
output_dir/splice/
output_dir/chunk_plan.json
```

**To resume a partial compute run**, keep the entire output directory
as-is (including `chunks_done/` and `splice/`).

#### Reclaiming disk space after a successful run

Once `alignable finalize` reports `complete` and you have what you
need, you can safely delete:

```bash
# Per-chunk raw artifacts (only needed for re-aggregation):
rm -r output_dir/splice

# If you only need the packaged zip:
rm -r output_dir/mappability.zarr     # ~27 GB; zip in sibling has the same data

# If you only need the unpackaged dir:
rm  output_dir.zarr.zip               # ~27 GB; dir has the same data
```

Keeping `config.yaml` in either case is recommended so you remember
how the store was built.

## Performance

The inner tiling loop uses a nanobind C++ extension (xoshiro256** PRNG,
batch read extraction, FASTQ formatting) achieving ~610K fragments/sec
with 1% error injection — a 33× speedup over pure Python. The scoring
loop uses a C++ extension backed by htslib for direct BAM parsing,
achieving a 250× speedup over the Python/pysam path.

| Configuration | chr22 (50 Mbp) | Est. full genome |
|---------------|----------------|------------------|
| STAR, 1 chunk, 24 threads | 26.6 min | ~16 hrs |
| STAR, 6 chunks parallel, shared memory | ~8 min | ~5 hrs |
| minimap2 splice:sr, 1 chunk, 24 threads | 42.5 min | ~26 hrs |

## CLI reference

### Commands

| Command | Description |
|---------|-------------|
| `alignable compute` | Full single-node run (init → compute → finalize) |
| `alignable init` | Initialize Zarr store and chunk plan only |
| `alignable compute-chunk` | Process a single chunk (for SLURM array jobs) |
| `alignable finalize` | Verify completeness and aggregate splice artifacts into `splice_blacklist.feather` |
| `alignable package` | Bundle the unpackaged `mappability.zarr/` into a standalone sibling `<output_dir>.zarr.zip` for transfer |
| `alignable query` | Query mappability for regions (BedGraph output) |
| `alignable summary` | Print store metadata and mean mappability |
| `alignable export-bigwig` | Export mappability as BigWig file |
| `alignable export-splice-bed` | Export `splice_blacklist.feather` as a gzipped BED file (opt-in; not produced by default) |

### Global options

| Flag | Description |
|------|-------------|
| `--config FILE` | YAML config file providing defaults for all parameters |
| `--version` | Show version |
| `--help` | Show help |

## License

See [LICENSE](LICENSE).
