Metadata-Version: 2.4
Name: arraysplitter
Version: 1.10.0
Classifier: Development Status :: 4 - Beta
Classifier: Environment :: Console
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: MacOS
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
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
License-File: LICENSE
Summary: De Novo Decomposition of Satellite DNA Arrays into Monomers within Telomere-to-Telomere Assemblies
Keywords: bioinformatics,satellite-dna,centromere,tandem-repeats,genomics
Author-email: Aleksey Komissarov <ad3002@gmail.com>
License: MIT
Requires-Python: >=3.8
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Homepage, https://github.com/aglabx/ArraySplitter
Project-URL: Issues, https://github.com/aglabx/ArraySplitter/issues
Project-URL: Repository, https://github.com/aglabx/ArraySplitter

# ArraySplitter: De Novo Decomposition of Satellite DNA Arrays

Decomposes satellite DNA arrays into monomers within telomere-to-telomere (T2T) assemblies. Ideal for analyzing centromeric and pericentromeric regions on monomeric level.

**Status:** Production ready. Successfully handles arrays from kilobase to megabase scale.

**Key Features:**
- De novo monomer identification without prior knowledge
- Autocorrelation-based period detection for robust periodicity analysis
- Automatic orientation to canonical form (A>T, C>G)
- Deterministic output sorted by genomic coordinates
- Multi-threaded processing

**Performance:** CHM13v2.0 assembly (~1300 alpha satellite arrays) processes in ~3.5 minutes (16 threads)

## Installation

```bash
pip install arraysplitter
```

Or build from source:
```bash
cd src/rust/arraysplitter
cargo build --release
```

## Quick Start

```bash
# Basic decomposition
arraysplitter -i arrays.fa -o output_prefix -t 16

# With predefined cut sequences
arraysplitter -i arrays.fa -o output_prefix -c ATG,CGCG -t 16

# Show version
arraysplitter --version
```

## Understanding Key Metrics

### `orientation` — Strand Direction
- **`fwd`** = sequence kept as-is (canonical form: more A's than T's)
- **`rev`** = sequence was reverse-complemented to canonical form

Why it matters: Satellite arrays on opposite strands appear as different sequences. Canonical orientation ensures consistent comparison across the genome.

### `ed_tmpl` — Edit Distance to Template (Consensus)
Measures how different each monomer is from the **consensus sequence** of all monomers in the array.

- **Low ed_tmpl (0-5)** = monomer is very similar to consensus → typical/canonical monomer
- **High ed_tmpl (>20)** = monomer is divergent → possibly a variant, mutation hotspot, or misalignment

Use case: Filter out divergent monomers, identify conserved vs variable positions.

### `ed_prev` / `ed_next` — Edit Distance to Neighbors
Measures how different each monomer is from its **adjacent monomers** (previous and next).

- **Low ed_prev/ed_next** = adjacent monomers are similar → homogeneous region
- **High ed_prev/ed_next** = sudden change → possible HOR boundary, structural variant, or sequence transition

Use case: Detect HOR structure boundaries, identify recombination breakpoints.

### `autocorr` — Autocorrelation Score
Measures periodicity strength (0.0 to 1.0):

- **High autocorr (>0.7)** = strong, regular periodicity → well-conserved tandem repeat
- **Medium autocorr (0.5-0.7)** = detectable but irregular periodicity
- **Low autocorr (<0.5)** = weak or no periodicity → degenerate or non-repetitive

### `cv` — Coefficient of Variation
Standard deviation of monomer lengths divided by mean length:

- **Low cv (<0.1)** = uniform monomer sizes → canonical repeat structure
- **High cv (>0.2)** = variable sizes → insertions/deletions, heterogeneous array

### `ed_per_bp` — Normalized Edit Distance
Edit distance divided by sequence length. Allows comparison across different monomer sizes:

- **0.00-0.02** = highly conserved
- **0.02-0.05** = moderately divergent
- **>0.05** = significantly divergent

## Output Files

All output is **deterministically sorted by chromosome and genomic position** (chr1 → chr22 → chrX → chrY → chrM).

| File | Description |
|------|-------------|
| `.decomposed.fasta` | Monomers with orientation info in headers |
| `.hors.tsv` | HOR-level decomposition (18 columns; includes level + parent_idx for intermediate HOR rows) |
| `.monomers.tsv` | Base-level monomers from recursive decomposition (18 columns; includes parent_level) |
| `.summary.tsv` | One-row-per-array summary with HOR and monomer statistics (28 columns) |
| `.templates.fasta` | Family templates built in the first pass (default mode; not written with `--templates` or `--no-phase-sync`) |
| `.lengths` | Fragment lengths for each array |

### Summary TSV Columns (`.summary.tsv`)

One row per array combining HOR-level and monomer-level statistics. Useful for overview analysis.

| Column | Description |
|--------|-------------|
| `array_id` | Array identifier (chr_start_end_len_period_type) |
| `array_length` | Total array length in bp |
| `orientation` | `fwd` or `rev` (reverse complemented to canonical) |
| `method` | Detection method used (`autocorr`, `classic`) |
| **HOR-level stats** | |
| `hor_period` | Detected HOR period in bp |
| `hor_autocorr` | Autocorrelation at HOR period |
| `hor_n_monomers` | Number of HOR-level monomers |
| `hor_mean_ed_tmpl` | Mean edit distance to HOR consensus |
| `hor_mean_ed_prev` | Mean edit distance between adjacent HORs |
| `hor_cv` | Coefficient of variation for HOR lengths |
| `hor_consensus` | Consensus sequence at HOR level |
| `hor_iupac` | IUPAC ambiguity codes (bases ≥20% frequency) |
| `hor_quality` | Per-position support (digit 0-9, 9=90-100%) |
| **Monomer-level stats** | |
| `mono_period` | Median base monomer period |
| `mono_autocorr` | Mean autocorrelation at monomer level |
| `mono_n_monomers` | Total number of base monomers |
| `mono_mean_ed_tmpl` | Mean edit distance to monomer consensus |
| `mono_mean_ed_prev` | Mean edit distance between adjacent monomers |
| `mono_cv` | Mean coefficient of variation |
| `mono_consensus` | Consensus sequence at monomer level |
| `mono_iupac` | IUPAC ambiguity codes |
| `mono_quality` | Per-position support |
| `cut_sequence` | Anchor k-mer used for splitting |
| `period_classes` | Base monomer size classes: `period:count,...` sorted by count desc (e.g., `171:85` = homogeneous, `171:85,2420:3` = heterogeneous) |
| `status` | `ok`, or `no_period` when no period was found (the whole array is one piece) |
| **Family phase** | (see [Family phase across arrays](#family-phase-across-arrays-reads)) |
| `family_id` | Template the array was cut at; `-` with `--no-phase-sync` or when no template matched |
| `phase_source` | `own_anchor` (`--no-phase-sync`), `template`, `template_hor` (HOR units snapped to the template), `template_no_cut`, `no_template_match` (no template matched: own anchor used), `no_template_of_period` (only templates of another length matched: own anchor used) |
| `template_score` | Fraction of the array's 11-mers found in its best template; `-` with `--no-phase-sync` |

### HORs TSV Columns (`.hors.tsv`)

Contains the primary decomposition into HOR (Higher Order Repeat) monomers. Multiple rows per array.

**Row types** (in order, per array):
1. `pred_array` - Array-level prediction/header row (level = 1)
2. `flank` - Terminal fragments <70% of period (level = 1)
3. `monomer` - Top-level HOR units (level = 1, sorted by idx)
4. `sub_hor` - Intermediate HOR nodes from recursive autocorrelation (level >= 2; only present when recursion finds further periodicity at autocorr ≥ 0.5)
5. `array` - Summary statistics row (level = 1)
6. `consensus` - Consensus sequence row (level = 1)

| Column | Description |
|--------|-------------|
| `array_id` | Array identifier (chr_start_end_len_period_type) |
| `type` | `pred_array`, `monomer`, `flank`, `sub_hor`, `array`, `consensus` |
| `idx` | Index within `(level, array_id)`; for level=1 rows this is the monomer position; for `sub_hor` rows the per-array counter at that level |
| `length` | Sequence length in bp |
| `source` | Detection method: `anchor`, `split_2x`, `split_3x`, `left_flank`, `right_flank`, `recursive_anchor`, `recursive_split` |
| `ed_tmpl` | Edit distance to consensus template |
| `ed_prev` | Edit distance to previous monomer |
| `ed_next` | Edit distance to next monomer |
| `period` | Detected repeat period in bp |
| `autocorr` | Autocorrelation value at detected period |
| `n_expected` | Expected count of monomers (array_len / period); for `sub_hor` rows this is the number of children |
| `ed_per_bp` | Normalized edit distance (ed / length) |
| `cv` | Coefficient of variation for lengths |
| `cut_sequence` | Anchor sequence used for splitting |
| `orientation` | `fwd` or `rev` (reverse complemented) |
| `sequence` | Actual DNA sequence (or `-` for pred_array/array rows) |
| `level` | 1 for top-level HOR units; 2, 3, ... for sub_hor rows by recursion depth |
| `parent_idx` | For `sub_hor` rows: idx within `level - 1` of the parent HOR. `-1` at level 1 (no parent above the array) |

### Monomers TSV Columns (`.monomers.tsv`)

Contains base-level monomers after recursive HOR decomposition. **Unified format** matching `.hors.tsv` plus `parent_idx`.

Each HOR is recursively decomposed until:
- No further periodicity detected (autocorrelation ≤ 0.5)
- Minimum length (5bp) reached

**Row types** (in order):
1. `pred_array` - Array-level summary row
2. `base_monomer` - Base-level monomers from recursive decomposition
3. `monomer` - Non-decomposable monomers (e.g., telomeres)

| Column | Description |
|--------|-------------|
| `array_id` | Array identifier |
| `type` | `pred_array`, `base_monomer`, `monomer` |
| `idx` | Global index within array (0-based) |
| `length` | Sequence length in bp |
| `source` | `recursive_anchor`, `recursive_split`, `base`, `recursive_flank` |
| `ed_tmpl` | Edit distance to submonomer consensus |
| `ed_prev` | Edit distance to previous base monomer |
| `ed_next` | Edit distance to next base monomer |
| `period` | Detected period at this level (0 if base) |
| `autocorr` | Autocorrelation value |
| `n_expected` | Always 1 for individual monomers |
| `ed_per_bp` | Normalized edit distance |
| `cv` | Coefficient of variation within parent group |
| `cut_sequence` | Inherited anchor sequence |
| `orientation` | Inherited from array (`fwd`/`rev`) |
| `parent_idx` | idx (within `parent_level`) of the deepest enclosing row in `.hors.tsv`. For leaves under a top-level HOR this is the top hor_idx; for leaves under a `sub_hor` this is that sub_hor's idx |
| `sequence` | Actual DNA sequence |
| `parent_level` | Level of the deepest enclosing HOR row: 1 if the leaf hangs directly off a top-level HOR; 2+ if it hangs off an intermediate `sub_hor` at that level. Walking `parent_idx` up `.hors.tsv` from this row until `level = 1` recovers the full ancestry chain |

### Example: α-satellite HOR Decomposition

For a typical α-satellite HOR (512bp → 3×171bp monomers):

**`.hors.tsv`** - 10 HOR monomers (~512bp each):
```
array_id                type        idx  length  period  ...
chr1_centromere         pred_array  10   5120    512     ...
chr1_centromere         monomer     0    512     512     ...
chr1_centromere         monomer     1    512     512     ...
...
chr1_centromere         array       10   5120    512     ...
chr1_centromere         consensus   10   512     512     ... [consensus seq]
```

**`.monomers.tsv`** - 30 base monomers (~171bp each):
```
array_id                type          idx  length  parent_idx  ...
chr1_centromere         pred_array    30   5120    -           ...
chr1_centromere         base_monomer  0    171     0           ...
chr1_centromere         base_monomer  1    171     0           ...
chr1_centromere         base_monomer  2    170     0           ...
chr1_centromere         base_monomer  3    171     1           ...
...
```

**`.summary.tsv`** - Single row with both levels:
```
array_id         length  hor_period  hor_n_monomers  mono_period  mono_n_monomers  period_classes  ...
chr1_centromere  5120    512         10              171          30               171:30          ...
```

## Algorithm

ArraySplitter employs an autocorrelation-based algorithm for detecting repeat periods and decomposing satellite DNA arrays.

### 1. Canonical Orientation

Arrays are oriented to canonical form:
- Primary rule: A > T (more A's than T's)
- Secondary rule: C > G (if A=T)
- Non-canonical arrays are reverse complemented
- With family phase (default), arrays matched to a template are oriented to the template strand instead

### 2. Period Detection via Autocorrelation

The algorithm computes sequence autocorrelation to detect periodicity:

```
autocorr(offset) = matches / comparisons
```

Where `matches` counts identical nucleotides at positions `i` and `i + offset`.

**Key innovations:**
- **Random expectation correction:** Subtracts expected random match rate based on nucleotide composition
- **Refined period search:** Harmonic folding peak detection to distinguish fundamental period vs harmonics
- **Confidence scoring:** Autocorrelation excess over random indicates detection confidence

### 3. Anchor Selection

For the detected period, finds optimal anchor (cut sequence) using:

1. **K-mer enumeration:** Evaluate candidate k-mers with k ∈ [3, min(11, period − 1)] from the sequence
2. **Position analysis:** For each k-mer, record all occurrence positions
3. **Scoring metrics:**
   - **Uniqueness:** Fraction of occurrences exactly `period` apart
   - **Regularity:** How evenly spaced the occurrences are
4. **Combined score:** `uniqueness × regularity`
5. **Deterministic selection:** K-mers sorted lexicographically for reproducible tie-breaking

### 4. Array Decomposition

Using the selected anchor:
1. Split array at all anchor occurrences
2. First fragment → left flank (if < 70% of period)
3. Middle fragments → monomers
4. Last fragment → right flank (if < 70% of period)
5. Apply heuristics for multiplet splitting (doublets, triplets, etc.)

### 5. Output Generation

Results are:
- Sorted by chromosome (natural order: 1, 2, ..., 22, X, Y, M)
- Within chromosome, sorted by start position
- Fully deterministic across runs

## Family phase across arrays (reads)

Without templates every array picks its own anchor k-mer: the best uniqueness × regularity score, ties broken
alphabetically. Two reads of the same locus can therefore be cut at different positions of the monomer (on a
simulated 3 × 178-bp HOR, 12 reads of one locus gave leaf boundaries at 20 different offsets; see
`src/rust/arraysplitter/tests/phase_sync.rs`). Where adjacent monomers differ (HORs), pieces cut at different phases are different sequences, and rotating them
afterwards does not help: rotation commutes with cutting only when neighbours are identical, so it produces
chimeras of two neighbours.

A **template** fixes the phase for a family: one monomer sequence whose position 0 is the cut point. Since v1.10.0
this is the default.

```bash
# Default: build templates from the input itself (first pass), then cut with them
arraysplitter -i reads.fa -o out                   # writes out.templates.fasta

# Templates you already have (plain A/C/G/T, one monomer per record)
arraysplitter -i reads.fa -o out --templates cen178.fa

# v1.9.x behaviour: every array at its own anchor, recursive sub_hor levels
arraysplitter -i reads.fa -o out --no-phase-sync
```

For each array:
1. The family and strand are chosen by shared 11-mers (`--template-min-score`, default 0.05), among templates
   within 0.8–1.25× of the array's own leaf period (the median leaf length of its own decomposition). A 532-bp
   trimer template shares more 11-mers with HOR reads than the 178-bp monomer template, but cutting 178-bp monomers
   on its lattice would make HOR units the leaves; an array that matches only templates of another length keeps its
   own anchor and is marked `no_template_of_period`. A matched array is oriented to the template strand (this
   replaces the A > T rule for matched arrays).
2. Every 11-mer shared with the template votes for the template position 0 of its monomer. Vote peaks are taken
   strongest first; a peak next to an accepted cut (within 1.5 template lengths) must be one template length from
   it (±15%), and a peak with no accepted cut that close starts a new stretch. The lattice phase is local, so an
   insertion or deletion in one monomer does not lose the cuts on the far side of it.
3. The pieces before the first cut and after the last cut are flanks (`template_left_flank` /
   `template_right_flank`). When the array's own period is a HOR (≥ 1.5 template lengths), its HOR boundaries are
   moved to the nearest cut, so HOR units also start at template position 0 (`phase_source = template_hor`).
   Which monomer opens the HOR is the same in every read of a locus: with m = HOR period / template length, the
   residue (mod m) of monomers that are on average most similar to the template marks the HOR start, and every HOR
   boundary moves to the nearest monomer of that residue.
4. Leaves (`.monomers.tsv`) are the template units themselves: every leaf runs from one position 0 to the next,
   including whole monomers inside the flank pieces. A unit ≥ 1.5 template lengths (one cut had too few votes) is
   kept whole and labelled `template_gap`.

Things to know:
- The template defines the leaf unit, so template-cut arrays have no recursive `sub_hor` rows.
- `cut_sequence` is the template's first 11 bp for every array of the family; `.decomposed.fasta` headers gain
  `family=` and `phase=`.
- Arrays that match no template are cut at their own anchor and marked `no_template_match`; the run log reports
  how many.
- What stays uncut: a stretch with no 11-mer shared with the template (a diverged variant) becomes one
  `template_gap` unit, as does the stretch after a deletion of more than half a monomer (a cut closer than half a
  template length to another is never kept) or a run of monomers more than 15% longer or shorter than the template.
  On the zebra finch panel `template_gap` units hold 6.4% of phased bp, 96% of it in stretches sharing < 5% of
  11-mers with the template; on *A. thaliana* reads 1.4%.
- Templates must be at least 22 bp, so short-period repeats (telomeric TTAGGG, microsatellites) get no family in
  the default mode and stay on their own anchor (`no_template_match`). On the zebra finch panel this was 99 of 429
  arrays, 75 of them period 6; the other 330 were phased into 30 families.
- The default first pass takes one medoid leaf per array, groups them by shared 11-mers (a candidate scoring < 0.25 to every
  family founds a new one) and writes them to `<out>.templates.fasta`. A new family takes its strand and position 0
  from its closest family (score ≥ 0.05), chosen among families of compatible length (0.8–1.25×) first, by the most
  common offset of shared 11-mers. A 138-bp deletion variant therefore cannot shift the cut of full 178-bp
  monomers. Position 0 of an automatic family is its seed array's own cut, so it is arbitrary but
  shared. To compare runs or datasets, cut them all with the same `--templates` file; rerunning with the written
  file reproduces the default output byte for byte. The first pass keeps each array's own decomposition, so the
  second pass does not recompute it; on the zebra finch panel the default run takes 1.27× the `--no-phase-sync`
  time (48.9 vs 38.5 s on 4 threads).
- If no array yields a usable unit, the run continues with own anchors, every array is marked `no_template_match`
  and the log says so.

## Methods

### `autocorr` (Default)

Uses autocorrelation for period detection. Best for:
- Regular tandem repeats
- Alpha satellite arrays
- HOR (Higher Order Repeat) structures

### `classic`

Uses frequency suffix tree approach. Better for:
- Irregular or degenerate repeats
- Very short arrays
- Arrays with high mutation rates

### `both`

Tries autocorrelation first, falls back to classic if autocorr fails.

## Command Line Options

```
arraysplitter --help

Options:
  -i, --input <FILE>       Input FASTA file
  -o, --output <PREFIX>    Output prefix
  -t, --threads <N>        Number of threads [default: all cores]
  -c, --cuts <SEQ,SEQ>     Predefined cut sequences (comma-separated)
  -d, --depth <N>          Max depth for cut search [default: 100]
  --method <METHOD>        Detection method: autocorr, classic, both [default: autocorr]
  --max-ed-len <N>         Max monomer length for edit distance [default: 10000]
  --max-period <N>         Max candidate period for autocorrelation scan in bp [default: 5000].
                           Larger values (e.g. 32000) surface very long HOR-level structure
                           that the default cap hides. Cost is O(L × max_period); on the zebra
                           finch panel raising to 32000 was ~33% slower wall-time and exposed
                           a different hor_period in 14 of 429 arrays. autocorr method only.
  --templates <FASTA>      Family templates: cut every matching array at template position 0
  --no-phase-sync          Own anchor per array, recursive sub_hor levels (v1.9.x behaviour)
  --phase-sync             Default; accepted for old command lines
  --template-min-score <F> Minimum 11-mer fraction to assign an array to a template [default: 0.05]
  --stats                  Print detailed statistics
  --top-outliers <N>       Number of outliers to show [default: 10]
  -V, --version            Print version
```

## Releases

- **v1.10.0** makes family phase the default. On 58,639 *A. thaliana* CEN178 reads every CEN178 family is cut at
  one phase, within 1 bp for 178- vs 179-bp variants (the four phases reported for the first draft are gone), and 76.6% of monomer occurrences fall in forms
  seen at least ½·c0 = 41 times (v1.9.1 own anchors: 62.7%). The default run takes 1.27× the `--no-phase-sync`
  time on the zebra finch panel.
  Family phase: a first pass builds family templates from the input, and every array
  (read) of a family is oriented to one strand and cut at one shared position of the monomer, so pieces of the same
  locus are comparable across reads. `--templates` supplies the templates instead; `--no-phase-sync` restores the
  v1.9.x cuts byte for byte (plus three new `.summary.tsv` columns: `family_id`, `phase_source`,
  `template_score`). Phased arrays have no `sub_hor` rows. The benchmark scripts and the paper notebook pass
  `--no-phase-sync` to binaries that support it, so they keep reproducing the published tables.
- **v1.9.1** keeps the piece before the first constellation cut as a flank child (v1.9.0 dropped it, so a few HOR
  copies were not tiled by their monomers), and adds `scripts/benchmark/check_tiling.py`, one-to-one boundary matching
  in the sweep-set scorer, and the corrected cat benchmark scripts (Table 5 / Figure 2 from one table, boundary check
  against `nhmmer`).
- **v1.9.0** adds a per-level regime router: at every level of the descent, arrays whose units are
  homogeneous keep the v1.8.2 single-anchor cut, and arrays whose monomers are too diverged for one exact
  anchor are cut with a circularly phased 7-mer constellation (`src/rust/arraysplitter/src/constellation.rs`).
  Benchmark: `scripts/benchmark/cat113/` (domestic cat 113-bp HORs).
- **v1.8.2** (commit `433722a`) is the single-anchor release and the baseline in that benchmark.

## Citation

If you use ArraySplitter in your research, please cite:
[Publication pending]

## Contact

For questions or support: ad3002@gmail.com

