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

pip install arraysplitter

Or build from source:

cd src/rust/arraysplitter
cargo build --release

Quick Start

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

Summary TSV Columns (.summary.tsv)

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

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)

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)

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.

# 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