From 51a9e2574d04b9f1a0b1497f757a13c6550b15b6 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Sun, 9 Aug 2026 16:19:50 +0200 Subject: [PATCH] feat: add stats and compare subcommands - stats: single streaming pass summarizing alignment amount, within-chain gaps, query strands, and chain fragmentation (aligned-bp N50/L50, per-target top-1/top-3 concentration, chains-to-50/90/95%, median/max fragmentation). --by-sequence adds per-target tables and top scored chains; reads stdin when --chain is omitted. - compare: canonical base-pair mapping agreement invariant to chain IDs, order, and harmless splitting, reported separately from unique coverage, multi-mapped bp, continuity, and gaps. Name interning plus a fixed 8 GiB memory budget checked before canonicalization; --by-sequence adds per-target agreement and top chains. - Bump to 0.0.11, add tool docs (assets/tools/stats.md, compare.md), README statistics section, and changelog entries. --- Cargo.lock | 2 +- Cargo.toml | 2 +- README.md | 9 +- assets/changelog/Changelog.md | 24 + assets/tools/compare.md | 155 ++++++ assets/tools/stats.md | 140 ++++++ src/cli/compare.rs | 916 ++++++++++++++++++++++++++++++++++ src/cli/filter.rs | 12 +- src/cli/merge.rs | 31 +- src/cli/mod.rs | 11 + src/cli/sort.rs | 13 +- src/cli/split.rs | 10 +- src/cli/stats.rs | 592 ++++++++++++++++++++++ src/io/index.rs | 2 +- src/io/reader.rs | 2 +- tests/chains.rs | 1 + 16 files changed, 1886 insertions(+), 36 deletions(-) create mode 100644 assets/tools/compare.md create mode 100644 assets/tools/stats.md create mode 100644 src/cli/compare.rs create mode 100644 src/cli/stats.rs diff --git a/Cargo.lock b/Cargo.lock index 7c1dd94..f23838b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -134,7 +134,7 @@ checksum = "f079e83a288787bcd14a6aea84cee5c87a67c5a3e660c30f557a3d24761b3527" [[package]] name = "chaintools" -version = "0.0.10" +version = "0.0.11" dependencies = [ "clap", "flate2", diff --git a/Cargo.toml b/Cargo.toml index 3c8bab2..e858f92 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "chaintools" -version = "0.0.10" +version = "0.0.11" edition = "2024" authors = ["alejandrogzi "] license = "MIT" diff --git a/README.md b/README.md index e701d17..f8bb6ee 100644 --- a/README.md +++ b/README.md @@ -10,8 +10,8 @@

- - Version Badge + + Version Badge Crates.io Version @@ -29,7 +29,7 @@ work with .chain files in Rust

- docs . + docs . usage . tools . chains @@ -60,8 +60,9 @@ Add this to your `Cargo.toml`: ```toml [dependencies] -chaintools = { version = "0.0.10", features = ["mmap", "gzip", "parallel"] } +chaintools = { version = "0.0.11", features = ["mmap", "gzip", "parallel"] } ``` + ## Benchmarks

diff --git a/assets/changelog/Changelog.md b/assets/changelog/Changelog.md index cb7269c..55d7eba 100644 --- a/assets/changelog/Changelog.md +++ b/assets/changelog/Changelog.md @@ -5,6 +5,30 @@ All notable changes to **chaintools** are documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.0.11] - 2026-08-09 + +### Added + +- **`stats` subcommand** — summarizes alignment amount, within-chain gaps, + query-strand counts, and chain fragmentation in one streaming pass. + Includes chain aligned-bp N50/L50, per-target top-1/top-3 chain + concentration and chains→50/90/95% counts, plus `--by-sequence` per-target + tables and top scored chains. Reads standard input when `--chain` is + omitted. + +- **`compare` subcommand** — compares two chain files by their canonical + base-pair mappings (invariant to chain IDs, order, and harmless splitting), + reported separately from unique coverage, multi-mapped bp, continuity, and + gaps. Per-target agreement is available with `--by-sequence`. Sequence names + are interned and a fixed memory budget is checked before canonicalization, + so oversized inputs fail cleanly instead of OOMing. + +### Documentation + +- Tool docs for the two new subcommands in `assets/tools/` + (`stats.md`, `compare.md`). +- README "Statistics and comparison" section with both commands. + ## [0.0.10] - 2026-08-08 ### Added diff --git a/assets/tools/compare.md b/assets/tools/compare.md new file mode 100644 index 0000000..eca4826 --- /dev/null +++ b/assets/tools/compare.md @@ -0,0 +1,155 @@ +

+

+ +

+ + +

+ chaintools +

+ + +

+ + Version Badge + + + Crates.io Version + + + GitHub License + + + Crates.io Total Downloads + +

+

+ + +# chaintools compare + +Compare two chain files by the exact base-pair relations they encode, separate +from how those relations are packaged into chains. + +## Usage + +```bash +chaintools compare --chain-a A.chain --chain-b B.chain +chaintools compare -a A.chain -b B.chain --by-sequence --top 3 +``` + +Both inputs must be paths; there is no stdin mode. + +## Options + +- `--chain-a ` / `-a`: first chain file. +- `--chain-b ` / `-b`: second chain file. +- `--by-sequence`: add the per-target agreement table and top-chain lists. +- `--top `: with `--by-sequence`, show the top `N` scored chains per target. + Default `3`. + +## What is compared + +Every aligned block is canonicalized into the base-to-base relation it +represents: a `(target, query, strand, diagonal)` key plus a target interval. +Within a key, overlapping and adjacent intervals are merged and duplicates are +removed. The result is independent of chain IDs, chain order, chain splitting, +and block splitting along the same mapping diagonal. + +This means a mapping split across two chains compares equal to the same mapping +in one chain — but the two files' *continuity* differs, and `compare` reports +both facts independently. Chain counts, header overlap, and approximate +coordinate matching are never used as equivalence tests. + +## Output + +- **MAPPING AGREEMENT** — unique pair bp in each file, shared bp, A-only and + B-only bp, `A_retained_by_B = shared / A`, `B_retained_by_A = shared / B`, + pair Dice and pair Jaccard. The directional values are not called + precision/recall: neither file is ground truth. +- **COVERAGE / AMBIGUITY** — unique target and query bp covered by the + mappings, and bp mapped more than once on each side. This distinguishes + genuinely broader coverage from redundant multi-mapping. +- **CONTINUITY** — chains, blocks, chain aligned-bp N50/L50, largest chain, + median blocks per chain, and chromosome fragmentation (median target top-1 + and top-3 share, median/max chains→90%), reusing the `stats` metrics. +- **GAPS** — target/query gap bp, gap-event counts, and largest gaps as A/B + rows with absolute and percent delta. Fewer gaps is not automatically + reported as better. +- **INTERPRETATION** — one or two lines drawn only from the measured + dimensions: mappings identical (possibly with a fragmentation difference), + or differing (Dice quoted, no winner declared); broader coverage with no + additional multi-mapping if that holds. There is no combined quality score. +- **BY TARGET** (with `--by-sequence`) — per-target continuity (chains, + top-1/top-3 share, N50, chains→90%) alongside per-target A-only/B-only/ + shared bp and Dice, so divergences can be localized. Followed by the top + `N` scored chains of each file per target. + +## Example: identical mappings, different fragmentation + +A holds one 24-bp chain; B holds the same two 12-bp mappings as two chains: + +```bash +$ chaintools compare -a frag_a.chain -b frag_b.chain +MAPPING AGREEMENT +A_unique_pair_bp 24 +B_unique_pair_bp 24 +shared_pair_bp 24 +A_only_pair_bp 0 +B_only_pair_bp 0 +A_retained_by_B 100.00% +B_retained_by_A 100.00% +pair_Dice 100.00% +pair_Jaccard 100.00% + +COVERAGE / AMBIGUITY +metric A B +unique_target_covered_bp 24 24 +unique_query_covered_bp 24 24 +target_bp_mapped_more_than_once 0 0 +query_bp_mapped_more_than_once 0 0 + +CONTINUITY +metric A B +chains 1 2 +blocks 2 2 +chain_aligned_bp_n50 24 12 +chain_aligned_bp_l50 1 1 +largest_chain_aligned_bp 24 12 +median_blocks_per_chain 2.0 1.0 +median_target_top1_share 100.00% 50.00% +median_target_top3_share 100.00% 100.00% +median_target_chains_to_90pct 1.0 2.0 +max_target_chains_to_90pct 1 2 + +GAPS +metric A B absolute_delta percent_delta +target_gap_bp 0 0 +0 0.00% +query_gap_bp 0 0 +0 0.00% +target_only_gap_events 0 0 +0 0.00% +query_only_gap_events 0 0 +0 0.00% +dual_gap_events 0 0 +0 0.00% +largest_target_gap 0 0 +0 0.00% +largest_query_gap 0 0 +0 0.00% + +INTERPRETATION +Canonical mappings are identical, but B is more fragmented. +``` + +100% mapping agreement and worse continuity in the same report: that is the +point of the tool. + +## Compatibility + +Before comparing, both files are fully validated (same checks as `stats`) and +sequence sizes are cross-checked. The same target or query name assigned +conflicting sizes in the two files is fatal. A sequence present in only one +file is not an error — that may be a real coverage difference. + +## Memory + +Canonical mappings are held in memory (one ~32-byte segment per block, both +files). The estimate is checked against a fixed 8 GiB ceiling before +canonicalization; past it the command fails with a clear error instead of +OOMing. Per-key external sorting is the planned upgrade path for larger inputs +and is not implemented. diff --git a/assets/tools/stats.md b/assets/tools/stats.md new file mode 100644 index 0000000..b1164df --- /dev/null +++ b/assets/tools/stats.md @@ -0,0 +1,140 @@ +

+

+ +

+ + +

+ chaintools +

+ + +

+ + Version Badge + + + Crates.io Version + + + GitHub License + + + Crates.io Total Downloads + +

+

+ + +# chaintools stats + +Summarize alignment amount, gaps, strands, and chain fragmentation in one +streaming pass over a chain file. + +## Usage + +```bash +chaintools stats --chain file.chain +chaintools stats --chain file.chain --by-sequence --top 3 +chaintools stats < file.chain +``` + +If `--chain` is omitted, chain data is read from standard input. `.gz` input +paths are supported when the `gzip` feature is enabled. + +## Options + +- `--chain `: input `.chain` file. Omit to read standard input. +- `--by-sequence`: add the per-target fragmentation table and top-chain lists. +- `--top `: with `--by-sequence`, show the top `N` scored chains per target. + Default `3`. + +## Output + +Output is tab-separated, grouped into sections: + +- **GENERAL** — chain and block counts, distinct target/query sequence counts. +- **ALIGNMENT** — `aligned_pair_bp` (sum of block sizes), `target_span_bp` + and `query_span_bp` (sum of header spans). `aligned_pair_bp` is never called + genome coverage. +- **CONTINUITY** — largest chain, chain aligned-bp N50/L50, median chain + aligned bp, median blocks per chain, and global chromosome fragmentation + (median/max target chains→90%). +- **GAPS** — within-chain gap bp on each side, gap-event types + (target-only `dt>0,dq==0`, query-only `dt==0,dq>0`, dual `dt>0,dq>0`), and + largest gaps. +- **STRAND** — plus/minus chain and aligned-bp counts (query strand). +- **SCORE** — sum of header scores and mean score per chain. +- **BY TARGET** (with `--by-sequence`) — one row per target: chain count, + total aligned bp, top-1 and top-3 aligned-bp share, chains required for + 50/90/95% of the target's aligned bp, and per-target chain N50. +- **TOP CHAINS** (with `--by-sequence`) — for each target, the top `N` chains + by score with aligned bp, spans, block count, alignment density + (`aligned_bp / target_span`), gaps, and largest gap. + +### Definitions + +- **top-1 share** = largest chain aligned bp / total chain aligned bp on the + target. **top-3 share** is the same with the largest 3 chains. These are + chain-concentration metrics, not unique genome coverage: one 100-Mb chain + scores 100%, two 50-Mb chains score 50%, even with identical base-pair + mappings. +- **chains→N%** = smallest `k` such that the `k` largest chains cover ≥ N% of + the target's aligned bp. +- **N50** uses chain aligned bp, not header span. +- Score alone is never used as a continuity metric; the top-chain table reports + density and spans alongside it. + +## Example + +```bash +$ chaintools stats --chain a.chain +GENERAL +chains 4 +blocks 7 +target_sequences 2 +query_sequences 3 + +ALIGNMENT +aligned_pair_bp 70 +target_span_bp 76 +query_span_bp 78 + +CONTINUITY +largest_chain_aligned_bp 40 +chain_aligned_bp_n50 40 +chain_aligned_bp_l50 1 +median_chain_aligned_bp 10.0 +median_blocks_per_chain 1.0 +median_target_top1_share 65.00% +median_target_top3_share 100.00% +median_target_chains_to_90pct 2.0 +max_target_chains_to_90pct 2 + +GAPS +target_gap_bp 6 +query_gap_bp 8 +target_only_gap_events 1 +query_only_gap_events 1 +dual_gap_events 1 +largest_target_gap 4 +largest_query_gap 5 + +STRAND +plus_chains 2 +minus_chains 2 +plus_aligned_bp 50 +minus_aligned_bp 20 + +SCORE +score_sum 200 +mean_score_per_chain 50.00 +``` + +## Validation + +Every chain is validated before it is counted: the target strand must be `+`, +spans must be non-empty and fit inside their sequence sizes, no block may be +empty, and block sizes plus gaps must account for both header spans exactly. +Conflicting sizes for the same sequence name anywhere in the input are fatal. +A malformed chain aborts with a `format error` naming the byte offset. diff --git a/src/cli/compare.rs b/src/cli/compare.rs new file mode 100644 index 0000000..82c099b --- /dev/null +++ b/src/cli/compare.rs @@ -0,0 +1,916 @@ +// Copyright (c) 2026 Alejandro Gonzales-Irribarren +// Distributed under the terms of the Apache License, Version 2.0. + +use std::cmp::Ordering; +use std::collections::{BTreeMap, BTreeSet, HashMap}; +use std::io::{BufRead, Write}; +use std::path::{Path, PathBuf}; + +use chaintools::{Strand, StreamingReader}; +use clap::Args; + +use super::CliError; +use super::stats::{self, ChainStats, TargetStats}; + +// ponytail: fixed in-memory ceiling; add per-key external sort if real inputs exceed 8 GiB. +const CANONICAL_MEMORY_BUDGET: u128 = 8 * 1024 * 1024 * 1024; +const ESTIMATED_SEGMENT_BYTES: u128 = 32; + +#[derive(Debug, Args)] +pub struct CompareArgs { + #[arg(short = 'a', long = "chain-a", value_name = "PATH")] + chain_a: PathBuf, + + #[arg(short = 'b', long = "chain-b", value_name = "PATH")] + chain_b: PathBuf, + + #[arg(long, help = "Add per-target agreement, fragmentation, and top chains")] + by_sequence: bool, + + #[arg( + long, + value_name = "N", + default_value_t = 3, + help = "Top scored chains to show per target with --by-sequence" + )] + top: usize, +} + +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)] +struct MappingKey { + target: u32, + query: u32, + strand: u8, + diagonal: i64, +} + +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)] +struct Segment { + key: MappingKey, + start: u32, + end: u32, +} + +#[derive(Default)] +struct NameInterner { + target_ids: HashMap, u32>, + query_ids: HashMap, u32>, + target_names: Vec>, + query_names: Vec>, +} + +impl NameInterner { + fn target(&mut self, name: &[u8]) -> Result { + intern_name(&mut self.target_ids, &mut self.target_names, name) + } + + fn query(&mut self, name: &[u8]) -> Result { + intern_name(&mut self.query_ids, &mut self.query_names, name) + } +} + +#[derive(Debug, Default, Clone, Copy)] +struct CoverageStats { + target_covered_bp: u64, + query_covered_bp: u64, + target_multi_bp: u64, + query_multi_bp: u64, +} + +#[derive(Debug, Default, Clone, Copy)] +struct TargetAgreement { + a_pair_bp: u64, + b_pair_bp: u64, + shared_pair_bp: u64, +} + +#[derive(Debug, Default)] +struct MappingAgreement { + a_pair_bp: u64, + b_pair_bp: u64, + shared_pair_bp: u64, + targets: Vec, +} + +#[derive(Debug, Clone, Copy)] +struct Event { + sequence: u32, + position: u32, + delta: i8, +} + +pub fn run( + args: CompareArgs, + _stdin: &mut R, + stdout: &mut W, + _stderr: &mut E, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, + E: Write, +{ + if args.top == 0 { + return Err(CliError::Message( + "--top must be greater than zero".to_owned(), + )); + } + super::ensure_inputs_exist( + &[("chain A", &args.chain_a), ("chain B", &args.chain_b)], + &[], + )?; + + let stats_a = analyze_path(&args.chain_a)?; + let stats_b = analyze_path(&args.chain_b)?; + ensure_compatible(&stats_a, &stats_b)?; + check_memory_budget(stats_a.blocks, stats_b.blocks)?; + + let mut names = NameInterner::default(); + let mappings_a = canonicalize_path(&args.chain_a, stats_a.blocks, &mut names)?; + let mappings_b = canonicalize_path(&args.chain_b, stats_b.blocks, &mut names)?; + let agreement = mapping_agreement(&mappings_a, &mappings_b, names.target_names.len()); + let coverage_a = coverage(&mappings_a)?; + let coverage_b = coverage(&mappings_b)?; + + write_comparison( + stdout, + &stats_a, + &stats_b, + &mappings_a, + &mappings_b, + &agreement, + coverage_a, + coverage_b, + &names, + args.by_sequence, + args.top, + )?; + stdout.flush()?; + super::log_summary( + "compare", + &[ + ("a_chains", stats_a.chains), + ("b_chains", stats_b.chains), + ("shared_pair_bp", agreement.shared_pair_bp), + ], + ); + Ok(()) +} + +fn analyze_path(path: &Path) -> Result { + stats::analyze(&mut StreamingReader::from_path(path)?) +} + +fn ensure_compatible(a: &ChainStats, b: &ChainStats) -> Result<(), CliError> { + check_sizes("target", &a.target_sizes, &b.target_sizes)?; + check_sizes("query", &a.query_sizes, &b.query_sizes) +} + +fn check_sizes( + side: &str, + a: &BTreeMap, u32>, + b: &BTreeMap, u32>, +) -> Result<(), CliError> { + for (name, size_a) in a { + if let Some(size_b) = b.get(name) + && size_a != size_b + { + return Err(CliError::Message(format!( + "conflicting {side} sizes for {}: A={size_a}, B={size_b}", + String::from_utf8_lossy(name) + ))); + } + } + Ok(()) +} + +fn check_memory_budget(a_blocks: u64, b_blocks: u64) -> Result<(), CliError> { + let estimated = (u128::from(a_blocks) + u128::from(b_blocks)) * ESTIMATED_SEGMENT_BYTES * 2; + if estimated > CANONICAL_MEMORY_BUDGET { + return Err(CliError::Message(format!( + "canonical mapping estimate ({:.2} GiB) exceeds compare's fixed 8 GiB ceiling; per-key external sort is the upgrade path and is not implemented", + estimated as f64 / (1024.0 * 1024.0 * 1024.0) + ))); + } + Ok(()) +} + +fn intern_name( + ids: &mut HashMap, u32>, + names: &mut Vec>, + name: &[u8], +) -> Result { + if let Some(id) = ids.get(name) { + return Ok(*id); + } + let id = u32::try_from(names.len()) + .map_err(|_| CliError::Message("more than 2^32 distinct sequence names".to_owned()))?; + let owned = name.to_vec(); + ids.insert(owned.clone(), id); + names.push(owned); + Ok(id) +} + +fn canonicalize_path( + path: &Path, + expected_blocks: u64, + names: &mut NameInterner, +) -> Result, CliError> { + canonicalize( + &mut StreamingReader::from_path(path)?, + expected_blocks, + names, + ) +} + +fn canonicalize( + reader: &mut StreamingReader, + expected_blocks: u64, + names: &mut NameInterner, +) -> Result, CliError> { + let capacity = usize::try_from(expected_blocks) + .map_err(|_| CliError::Message("block count exceeds addressable memory".to_owned()))?; + let mut segments = Vec::new(); + segments.try_reserve_exact(capacity).map_err(|_| { + CliError::Message("unable to reserve memory for canonical mapping segments".to_owned()) + })?; + + while let Some(header) = reader.next_header()? { + let offset = header.offset; + let blocks = reader.read_blocks(offset)?; + let chain = header.into_chain(blocks); + stats::validate_chain(&chain, offset)?; + let target = names.target(&chain.reference_name)?; + let query = names.query(&chain.query_name)?; + let strand = match chain.query_strand { + Strand::Plus => 0, + Strand::Minus => 1, + }; + let mut target_position = u64::from(chain.reference_start); + let mut query_position = u64::from(chain.query_start); + let last = chain.blocks.len().saturating_sub(1); + + for (index, block) in chain.blocks.iter().enumerate() { + let size = u64::from(block.size); + let target_end = target_position.checked_add(size).ok_or_else(|| { + CliError::Message("target block coordinate overflows u64".to_owned()) + })?; + let query_end = query_position.checked_add(size).ok_or_else(|| { + CliError::Message("query block coordinate overflows u64".to_owned()) + })?; + let diagonal = match chain.query_strand { + Strand::Plus => query_position as i64 - target_position as i64, + Strand::Minus => { + let query_forward_end = u64::from(chain.query_size) + .checked_sub(query_position) + .ok_or_else(|| { + CliError::Message( + "query block coordinate exceeds its sequence size".to_owned(), + ) + })?; + i64::try_from(target_position + query_forward_end).map_err(|_| { + CliError::Message("minus-strand mapping diagonal overflows i64".to_owned()) + })? + } + }; + segments.push(Segment { + key: MappingKey { + target, + query, + strand, + diagonal, + }, + start: u32::try_from(target_position) + .map_err(|_| CliError::Message("target block start exceeds u32".to_owned()))?, + end: u32::try_from(target_end) + .map_err(|_| CliError::Message("target block end exceeds u32".to_owned()))?, + }); + target_position = target_end; + query_position = query_end; + if index < last { + target_position += u64::from(block.gap_reference); + query_position += u64::from(block.gap_query); + } + } + } + + segments.sort_unstable(); + let mut kept = 0usize; + for index in 0..segments.len() { + let segment = segments[index]; + if kept > 0 + && segments[kept - 1].key == segment.key + && segment.start <= segments[kept - 1].end + { + segments[kept - 1].end = segments[kept - 1].end.max(segment.end); + } else { + segments[kept] = segment; + kept += 1; + } + } + segments.truncate(kept); + Ok(segments) +} + +fn mapping_agreement(a: &[Segment], b: &[Segment], target_count: usize) -> MappingAgreement { + let mut result = MappingAgreement { + targets: vec![TargetAgreement::default(); target_count], + ..MappingAgreement::default() + }; + for segment in a { + let length = u64::from(segment.end - segment.start); + result.a_pair_bp += length; + result.targets[segment.key.target as usize].a_pair_bp += length; + } + for segment in b { + let length = u64::from(segment.end - segment.start); + result.b_pair_bp += length; + result.targets[segment.key.target as usize].b_pair_bp += length; + } + + let (mut i, mut j) = (0usize, 0usize); + while i < a.len() && j < b.len() { + match a[i].key.cmp(&b[j].key) { + Ordering::Less => i += 1, + Ordering::Greater => j += 1, + Ordering::Equal if a[i].end <= b[j].start => i += 1, + Ordering::Equal if b[j].end <= a[i].start => j += 1, + Ordering::Equal => { + let shared = u64::from(a[i].end.min(b[j].end) - a[i].start.max(b[j].start)); + result.shared_pair_bp += shared; + result.targets[a[i].key.target as usize].shared_pair_bp += shared; + match a[i].end.cmp(&b[j].end) { + Ordering::Less => i += 1, + Ordering::Greater => j += 1, + Ordering::Equal => { + i += 1; + j += 1; + } + } + } + } + } + result +} + +fn coverage(segments: &[Segment]) -> Result { + let (target_covered_bp, target_multi_bp) = projected_coverage(segments, false)?; + let (query_covered_bp, query_multi_bp) = projected_coverage(segments, true)?; + Ok(CoverageStats { + target_covered_bp, + query_covered_bp, + target_multi_bp, + query_multi_bp, + }) +} + +fn projected_coverage(segments: &[Segment], query_side: bool) -> Result<(u64, u64), CliError> { + let capacity = segments + .len() + .checked_mul(2) + .ok_or_else(|| CliError::Message("coverage event count overflows usize".to_owned()))?; + let mut events = Vec::new(); + events + .try_reserve_exact(capacity) + .map_err(|_| CliError::Message("unable to reserve memory for coverage sweep".to_owned()))?; + for segment in segments { + let (sequence, start, end) = if query_side { + let (start, end) = query_interval(*segment)?; + (segment.key.query, start, end) + } else { + (segment.key.target, segment.start, segment.end) + }; + events.push(Event { + sequence, + position: start, + delta: 1, + }); + events.push(Event { + sequence, + position: end, + delta: -1, + }); + } + events.sort_unstable_by_key(|event| (event.sequence, event.position)); + + let mut covered = 0u64; + let mut multi = 0u64; + let mut current_sequence = None; + let mut previous_position = 0u32; + let mut depth = 0i64; + let mut index = 0usize; + while index < events.len() { + let sequence = events[index].sequence; + let position = events[index].position; + if current_sequence == Some(sequence) { + let length = u64::from(position - previous_position); + if depth > 0 { + covered += length; + } + if depth > 1 { + multi += length; + } + } else { + if current_sequence.is_some() && depth != 0 { + return Err(CliError::Message( + "unbalanced coverage intervals during sweep".to_owned(), + )); + } + current_sequence = Some(sequence); + depth = 0; + } + + let mut delta = 0i64; + while index < events.len() + && events[index].sequence == sequence + && events[index].position == position + { + delta += i64::from(events[index].delta); + index += 1; + } + depth += delta; + if depth < 0 { + return Err(CliError::Message( + "negative coverage depth during sweep".to_owned(), + )); + } + previous_position = position; + } + if depth != 0 { + return Err(CliError::Message( + "unbalanced coverage intervals during sweep".to_owned(), + )); + } + Ok((covered, multi)) +} + +fn query_interval(segment: Segment) -> Result<(u32, u32), CliError> { + let (start, end) = if segment.key.strand == 0 { + ( + i64::from(segment.start) + segment.key.diagonal, + i64::from(segment.end) + segment.key.diagonal, + ) + } else { + ( + segment.key.diagonal - i64::from(segment.end), + segment.key.diagonal - i64::from(segment.start), + ) + }; + if start < 0 || end <= start || end > i64::from(u32::MAX) { + return Err(CliError::Message( + "canonical mapping produced an invalid query interval".to_owned(), + )); + } + Ok((start as u32, end as u32)) +} + +#[allow(clippy::too_many_arguments)] +fn write_comparison( + writer: &mut W, + a: &ChainStats, + b: &ChainStats, + mappings_a: &[Segment], + mappings_b: &[Segment], + agreement: &MappingAgreement, + coverage_a: CoverageStats, + coverage_b: CoverageStats, + names: &NameInterner, + by_sequence: bool, + top: usize, +) -> Result<(), CliError> { + let a_only = agreement.a_pair_bp - agreement.shared_pair_bp; + let b_only = agreement.b_pair_bp - agreement.shared_pair_bp; + + writeln!(writer, "MAPPING AGREEMENT")?; + writeln!(writer, "A_unique_pair_bp\t{}", agreement.a_pair_bp)?; + writeln!(writer, "B_unique_pair_bp\t{}", agreement.b_pair_bp)?; + writeln!(writer, "shared_pair_bp\t{}", agreement.shared_pair_bp)?; + writeln!(writer, "A_only_pair_bp\t{a_only}")?; + writeln!(writer, "B_only_pair_bp\t{b_only}")?; + writeln!( + writer, + "A_retained_by_B\t{}", + percent_ratio(agreement.shared_pair_bp, agreement.a_pair_bp) + )?; + writeln!( + writer, + "B_retained_by_A\t{}", + percent_ratio(agreement.shared_pair_bp, agreement.b_pair_bp) + )?; + writeln!( + writer, + "pair_Dice\t{}", + percent_ratio_u128( + u128::from(agreement.shared_pair_bp) * 2, + u128::from(agreement.a_pair_bp) + u128::from(agreement.b_pair_bp) + ) + )?; + writeln!( + writer, + "pair_Jaccard\t{}", + percent_ratio_u128( + u128::from(agreement.shared_pair_bp), + u128::from(agreement.a_pair_bp) + u128::from(agreement.b_pair_bp) + - u128::from(agreement.shared_pair_bp) + ) + )?; + + writeln!(writer, "\nCOVERAGE / AMBIGUITY")?; + writeln!(writer, "metric\tA\tB")?; + writeln!( + writer, + "unique_target_covered_bp\t{}\t{}", + coverage_a.target_covered_bp, coverage_b.target_covered_bp + )?; + writeln!( + writer, + "unique_query_covered_bp\t{}\t{}", + coverage_a.query_covered_bp, coverage_b.query_covered_bp + )?; + writeln!( + writer, + "target_bp_mapped_more_than_once\t{}\t{}", + coverage_a.target_multi_bp, coverage_b.target_multi_bp + )?; + writeln!( + writer, + "query_bp_mapped_more_than_once\t{}\t{}", + coverage_a.query_multi_bp, coverage_b.query_multi_bp + )?; + + writeln!(writer, "\nCONTINUITY")?; + writeln!(writer, "metric\tA\tB")?; + writeln!(writer, "chains\t{}\t{}", a.chains, b.chains)?; + writeln!(writer, "blocks\t{}\t{}", a.blocks, b.blocks)?; + writeln!( + writer, + "chain_aligned_bp_n50\t{}\t{}", + a.chain_aligned_n50, b.chain_aligned_n50 + )?; + writeln!( + writer, + "chain_aligned_bp_l50\t{}\t{}", + a.chain_aligned_l50, b.chain_aligned_l50 + )?; + writeln!( + writer, + "largest_chain_aligned_bp\t{}\t{}", + a.largest_chain_aligned_bp, b.largest_chain_aligned_bp + )?; + writeln!( + writer, + "median_blocks_per_chain\t{:.1}\t{:.1}", + a.median_blocks_per_chain, b.median_blocks_per_chain + )?; + writeln!( + writer, + "median_target_top1_share\t{}\t{}", + percent(a.median_target_top1_share), + percent(b.median_target_top1_share) + )?; + writeln!( + writer, + "median_target_top3_share\t{}\t{}", + percent(a.median_target_top3_share), + percent(b.median_target_top3_share) + )?; + writeln!( + writer, + "median_target_chains_to_90pct\t{:.1}\t{:.1}", + a.median_target_chains_90, b.median_target_chains_90 + )?; + writeln!( + writer, + "max_target_chains_to_90pct\t{}\t{}", + a.max_target_chains_90, b.max_target_chains_90 + )?; + + writeln!(writer, "\nGAPS")?; + writeln!(writer, "metric\tA\tB\tabsolute_delta\tpercent_delta")?; + write_delta_row(writer, "target_gap_bp", a.target_gap_bp, b.target_gap_bp)?; + write_delta_row(writer, "query_gap_bp", a.query_gap_bp, b.query_gap_bp)?; + write_delta_row( + writer, + "target_only_gap_events", + a.target_only_gap_events, + b.target_only_gap_events, + )?; + write_delta_row( + writer, + "query_only_gap_events", + a.query_only_gap_events, + b.query_only_gap_events, + )?; + write_delta_row( + writer, + "dual_gap_events", + a.dual_gap_events, + b.dual_gap_events, + )?; + write_delta_row( + writer, + "largest_target_gap", + a.largest_target_gap, + b.largest_target_gap, + )?; + write_delta_row( + writer, + "largest_query_gap", + a.largest_query_gap, + b.largest_query_gap, + )?; + + write_interpretation(writer, a, b, agreement, coverage_a, coverage_b)?; + if by_sequence { + write_by_sequence(writer, a, b, agreement, names, top)?; + } + + debug_assert_eq!( + mappings_a.iter().map(segment_len).sum::(), + agreement.a_pair_bp + ); + debug_assert_eq!( + mappings_b.iter().map(segment_len).sum::(), + agreement.b_pair_bp + ); + Ok(()) +} + +fn write_delta_row(writer: &mut W, label: &str, a: u64, b: u64) -> Result<(), CliError> { + let delta = i128::from(b) - i128::from(a); + let percent_delta = if a == 0 { + if b == 0 { + "0.00%".to_owned() + } else { + "n/a".to_owned() + } + } else { + format!("{:.2}%", delta as f64 * 100.0 / a as f64) + }; + writeln!(writer, "{label}\t{a}\t{b}\t{delta:+}\t{percent_delta}")?; + Ok(()) +} + +fn write_interpretation( + writer: &mut W, + a: &ChainStats, + b: &ChainStats, + agreement: &MappingAgreement, + coverage_a: CoverageStats, + coverage_b: CoverageStats, +) -> Result<(), CliError> { + writeln!(writer, "\nINTERPRETATION")?; + let identical = agreement.a_pair_bp == agreement.shared_pair_bp + && agreement.b_pair_bp == agreement.shared_pair_bp; + if identical && more_fragmented(b, a) { + writeln!( + writer, + "Canonical mappings are identical, but B is more fragmented." + )?; + } else if identical && more_fragmented(a, b) { + writeln!( + writer, + "Canonical mappings are identical, but A is more fragmented." + )?; + } else if identical { + writeln!(writer, "Canonical base-pair mappings are identical.")?; + } else { + writeln!( + writer, + "Canonical mappings differ (pair Dice {}); continuity alone cannot identify a winner.", + percent_ratio_u128( + u128::from(agreement.shared_pair_bp) * 2, + u128::from(agreement.a_pair_bp) + u128::from(agreement.b_pair_bp) + ) + )?; + } + + if broader_without_more_ambiguity(coverage_b, coverage_a) { + writeln!( + writer, + "B has broader target/query coverage without more multi-mapped target/query bp." + )?; + } else if broader_without_more_ambiguity(coverage_a, coverage_b) { + writeln!( + writer, + "A has broader target/query coverage without more multi-mapped target/query bp." + )?; + } + Ok(()) +} + +fn more_fragmented(candidate: &ChainStats, other: &ChainStats) -> bool { + candidate.chain_aligned_n50 < other.chain_aligned_n50 + && candidate.median_target_top1_share < other.median_target_top1_share + && candidate.median_target_chains_90 > other.median_target_chains_90 +} + +fn broader_without_more_ambiguity(candidate: CoverageStats, other: CoverageStats) -> bool { + candidate.target_covered_bp >= other.target_covered_bp + && candidate.query_covered_bp >= other.query_covered_bp + && (candidate.target_covered_bp > other.target_covered_bp + || candidate.query_covered_bp > other.query_covered_bp) + && candidate.target_multi_bp <= other.target_multi_bp + && candidate.query_multi_bp <= other.query_multi_bp +} + +fn write_by_sequence( + writer: &mut W, + a: &ChainStats, + b: &ChainStats, + agreement: &MappingAgreement, + names: &NameInterner, + top: usize, +) -> Result<(), CliError> { + let target_names = a + .targets + .keys() + .chain(b.targets.keys()) + .cloned() + .collect::>(); + writeln!(writer, "\nBY TARGET")?; + writeln!( + writer, + "target\tA_chains\tB_chains\tA_top1\tB_top1\tA_top3\tB_top3\tA_N50\tB_N50\tA_chains_to_90pct\tB_chains_to_90pct\tA_only_pair_bp\tB_only_pair_bp\tshared_pair_bp\tpair_Dice" + )?; + for name in &target_names { + let target_a = a.targets.get(name); + let target_b = b.targets.get(name); + let pair = names + .target_ids + .get(name.as_slice()) + .and_then(|id| agreement.targets.get(*id as usize)) + .copied() + .unwrap_or_default(); + writeln!( + writer, + "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}", + String::from_utf8_lossy(name), + target_a.map_or(0, TargetStats::chain_count), + target_b.map_or(0, TargetStats::chain_count), + percent(target_a.map_or(0.0, |target| target.top1_share)), + percent(target_b.map_or(0.0, |target| target.top1_share)), + percent(target_a.map_or(0.0, |target| target.top3_share)), + percent(target_b.map_or(0.0, |target| target.top3_share)), + target_a.map_or(0, |target| target.chain_n50), + target_b.map_or(0, |target| target.chain_n50), + target_a.map_or(0, |target| target.chains_90), + target_b.map_or(0, |target| target.chains_90), + pair.a_pair_bp - pair.shared_pair_bp, + pair.b_pair_bp - pair.shared_pair_bp, + pair.shared_pair_bp, + percent_ratio_u128( + u128::from(pair.shared_pair_bp) * 2, + u128::from(pair.a_pair_bp) + u128::from(pair.b_pair_bp) + ) + )?; + } + + for name in &target_names { + if let Some(target) = a.targets.get(name) { + writeln!(writer, "\nTOP CHAINS A\t{}", String::from_utf8_lossy(name))?; + stats::write_top_chains(writer, &target.chain_summaries, top)?; + } + if let Some(target) = b.targets.get(name) { + writeln!(writer, "\nTOP CHAINS B\t{}", String::from_utf8_lossy(name))?; + stats::write_top_chains(writer, &target.chain_summaries, top)?; + } + } + Ok(()) +} + +fn segment_len(segment: &Segment) -> u64 { + u64::from(segment.end - segment.start) +} + +fn percent(value: f64) -> String { + format!("{:.2}%", value * 100.0) +} + +fn percent_ratio(numerator: u64, denominator: u64) -> String { + percent_ratio_u128(u128::from(numerator), u128::from(denominator)) +} + +fn percent_ratio_u128(numerator: u128, denominator: u128) -> String { + if denominator == 0 { + "n/a".to_owned() + } else { + format!("{:.2}%", numerator as f64 * 100.0 / denominator as f64) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use std::io::Cursor; + + fn stats_text(text: &str) -> ChainStats { + stats::analyze(&mut StreamingReader::new(Cursor::new(text.as_bytes()))) + .expect("valid stats fixture") + } + + fn canonical_text(text: &str, names: &mut NameInterner) -> Vec { + canonicalize( + &mut StreamingReader::new(Cursor::new(text.as_bytes())), + 16, + names, + ) + .expect("valid canonical fixture") + } + + #[test] + fn split_chains_keep_exact_mapping_but_expose_fragmentation() { + let a = "chain 100 chr1 100 + 0 24 qry1 100 + 0 24 7\n12 0 0\n12\n\n"; + let b = "chain 50 chr1 100 + 0 12 qry1 100 + 0 12 99\n12\n\n\ + chain 50 chr1 100 + 12 24 qry1 100 + 12 24 42\n12\n\n"; + let stats_a = stats_text(a); + let stats_b = stats_text(b); + let mut names = NameInterner::default(); + let mappings_a = canonical_text(a, &mut names); + let mappings_b = canonical_text(b, &mut names); + let agreement = mapping_agreement(&mappings_a, &mappings_b, names.target_names.len()); + let target_a = stats_a.targets.get(b"chr1".as_slice()).expect("A target"); + let target_b = stats_b.targets.get(b"chr1".as_slice()).expect("B target"); + + assert_eq!(agreement.a_pair_bp, 24); + assert_eq!(agreement.b_pair_bp, 24); + assert_eq!(agreement.shared_pair_bp, 24); + assert_eq!(target_a.top1_share, 1.0); + assert_eq!(target_b.top1_share, 0.5); + assert_eq!((target_a.chains_90, target_b.chains_90), (1, 2)); + assert!(more_fragmented(&stats_b, &stats_a)); + } + + #[test] + fn canonicalization_ignores_chain_ids_order_splits_and_duplicates() { + let a = "chain 1 chr1 100 + 0 10 q1 100 + 0 10 1\n10\n\n\ + chain 1 chr1 100 + 20 30 q1 100 + 20 30 2\n10\n\n"; + let b = "chain 9 chr1 100 + 20 30 q1 100 + 20 30 77\n10\n\n\ + chain 9 chr1 100 + 0 10 q1 100 + 0 10 88\n5 0 0\n5\n\n\ + chain 9 chr1 100 + 0 10 q1 100 + 0 10 99\n10\n\n"; + let mut names = NameInterner::default(); + let mappings_a = canonical_text(a, &mut names); + let mappings_b = canonical_text(b, &mut names); + + assert_eq!(mappings_a, mappings_b); + assert_eq!(mappings_a.iter().map(segment_len).sum::(), 20); + } + + #[test] + fn intersection_distinguishes_query_mapping_and_partial_overlap() { + let duplicated = "chain 1 chr1 100 + 0 10 q1 100 + 0 10 1\n10\n\n\ + chain 1 chr1 100 + 0 10 q1 100 + 0 10 2\n10\n\n"; + let shifted = "chain 1 chr1 100 + 5 15 q1 100 + 5 15 3\n10\n\n"; + let different_query = "chain 1 chr1 100 + 0 10 q1 100 + 1 11 4\n10\n\n"; + let mut names = NameInterner::default(); + let a = canonical_text(duplicated, &mut names); + let b = canonical_text(shifted, &mut names); + let c = canonical_text(different_query, &mut names); + let overlap = mapping_agreement(&a, &b, names.target_names.len()); + let mismatch = mapping_agreement(&a, &c, names.target_names.len()); + + assert_eq!( + (overlap.a_pair_bp, overlap.b_pair_bp, overlap.shared_pair_bp), + (10, 10, 5) + ); + assert_eq!(mismatch.shared_pair_bp, 0); + } + + #[test] + fn negative_strand_split_is_exact() { + let a = "chain 1 chr1 100 + 0 10 q1 100 - 20 30 1\n10\n\n"; + let b = "chain 1 chr1 100 + 5 10 q1 100 - 25 30 3\n5\n\n\ + chain 1 chr1 100 + 0 5 q1 100 - 20 25 2\n5\n\n"; + let mut names = NameInterner::default(); + let mappings_a = canonical_text(a, &mut names); + let mappings_b = canonical_text(b, &mut names); + let agreement = mapping_agreement(&mappings_a, &mappings_b, names.target_names.len()); + + assert_eq!(mappings_a, mappings_b); + assert_eq!(agreement.shared_pair_bp, 10); + } + + #[test] + fn compatibility_rejects_conflicting_sizes() { + let a = stats_text("chain 1 chr1 100 + 0 10 q1 100 + 0 10 1\n10\n\n"); + let b = stats_text("chain 1 chr1 101 + 0 10 q1 100 + 0 10 2\n10\n\n"); + let error = ensure_compatible(&a, &b).expect_err("size conflict"); + + assert!(error.to_string().contains("conflicting target sizes")); + } + + #[test] + fn coverage_sweep_counts_target_and_query_multimapping() { + let text = "chain 1 chr1 100 + 0 10 q1 100 + 0 10 1\n10\n\n\ + chain 1 chr1 100 + 0 10 q2 100 + 0 10 2\n10\n\n\ + chain 1 chr2 100 + 0 10 q1 100 + 0 10 3\n10\n\n"; + let mut names = NameInterner::default(); + let mappings = canonical_text(text, &mut names); + let coverage = coverage(&mappings).expect("coverage"); + + assert_eq!(coverage.target_covered_bp, 20); + assert_eq!(coverage.target_multi_bp, 10); + assert_eq!(coverage.query_covered_bp, 20); + assert_eq!(coverage.query_multi_bp, 10); + } +} diff --git a/src/cli/filter.rs b/src/cli/filter.rs index 775d027..c0c6a68 100644 --- a/src/cli/filter.rs +++ b/src/cli/filter.rs @@ -444,12 +444,12 @@ fn validate_output_args(args: &FilterArgs) -> Result<(), CliError> { /// Validates that output path differs from input paths. fn validate_output_path(args: &FilterArgs) -> Result<(), CliError> { - if let Some(output) = &args.out_chain { - if args.chains.iter().any(|input| input == output) { - return Err(CliError::Message( - "--out-chain must not be the same path as an input --chain".to_owned(), - )); - } + if let Some(output) = &args.out_chain + && args.chains.iter().any(|input| input == output) + { + return Err(CliError::Message( + "--out-chain must not be the same path as an input --chain".to_owned(), + )); } Ok(()) } diff --git a/src/cli/merge.rs b/src/cli/merge.rs index 97e6ab6..d012164 100644 --- a/src/cli/merge.rs +++ b/src/cli/merge.rs @@ -5,16 +5,16 @@ use std::fs::File; use std::io::{BufRead, BufReader, BufWriter, Read, Write}; use std::path::{Path, PathBuf}; -use chaintools::{io::storage::is_gz_path, StreamingReader}; +use chaintools::{StreamingReader, io::storage::is_gz_path}; use clap::{Args, ValueEnum}; #[cfg(feature = "gzip")] -use flate2::{read::MultiGzDecoder, write::GzEncoder, Compression}; +use flate2::{Compression, read::MultiGzDecoder, write::GzEncoder}; +use super::CliError; use super::sort_core::{ - emit_sorted_chains, write_metadata_lines, SortAccumulator, SortCriterion, - OUTPUT_BUFFER_CAPACITY, + OUTPUT_BUFFER_CAPACITY, SortAccumulator, SortCriterion, emit_sorted_chains, + write_metadata_lines, }; -use super::CliError; const BYTES_PER_GB: f64 = 1_000_000_000.0; const DEFAULT_MAX_GB: f64 = 16.0; @@ -793,9 +793,10 @@ mod tests { input.arg(), ]); - assert!(err - .to_string() - .contains("--out-chain must not be the same path as input chain")); + assert!( + err.to_string() + .contains("--out-chain must not be the same path as input chain") + ); } #[test] @@ -1042,9 +1043,10 @@ mod tests { arg("--gzip"), ]); - assert!(err - .to_string() - .contains("--gzip requires chaintools to be built with the `gzip` feature")); + assert!( + err.to_string() + .contains("--gzip requires chaintools to be built with the `gzip` feature") + ); } #[test] @@ -1064,8 +1066,9 @@ mod tests { arg("0"), ]); - assert!(err - .to_string() - .contains("--max-gb must be greater than zero")); + assert!( + err.to_string() + .contains("--max-gb must be greater than zero") + ); } } diff --git a/src/cli/mod.rs b/src/cli/mod.rs index 716d2fd..4e4aa90 100644 --- a/src/cli/mod.rs +++ b/src/cli/mod.rs @@ -3,6 +3,7 @@ pub mod anti_repeat; pub mod bed; +pub mod compare; pub mod convert; pub mod filter; pub mod liftover; @@ -11,6 +12,7 @@ pub mod score; pub mod sort; mod sort_core; pub mod split; +pub mod stats; pub mod swap; use std::fmt; @@ -56,6 +58,7 @@ pub struct Cli { /// Subcommands available in the chaintools CLI. #[derive(Debug, Subcommand)] +#[allow(clippy::large_enum_variant)] enum Command { #[command( name = "antirepeat", @@ -65,6 +68,8 @@ enum Command { AntiRepeat(anti_repeat::AntiRepeatArgs), #[command(about = "Convert chain files to BED")] Bed(bed::BedArgs), + #[command(about = "Compare exact mappings, coverage, ambiguity, and continuity")] + Compare(compare::CompareArgs), #[command(about = "Convert chain files to VCF, MAF, or BAM")] Convert(convert::ConvertArgs), #[command(about = "Filter chain files")] @@ -79,6 +84,8 @@ enum Command { Split(split::SplitArgs), #[command(about = "Sort chain files")] Sort(sort::SortArgs), + #[command(about = "Summarize alignment, gap, and continuity statistics")] + Stats(stats::StatsArgs), #[command(about = "Swap the target and query sides of chains (UCSC chainSwap)")] Swap(swap::SwapArgs), } @@ -88,6 +95,7 @@ impl std::fmt::Display for Command { match self { Command::AntiRepeat(_) => f.write_str("antirepeat"), Command::Bed(_) => f.write_str("bed"), + Command::Compare(_) => f.write_str("compare"), Command::Convert(_) => f.write_str("convert"), Command::Filter(_) => f.write_str("filter"), Command::Liftover(_) => f.write_str("liftover"), @@ -95,6 +103,7 @@ impl std::fmt::Display for Command { Command::Score(_) => f.write_str("score"), Command::Split(_) => f.write_str("split"), Command::Sort(_) => f.write_str("sort"), + Command::Stats(_) => f.write_str("stats"), Command::Swap(_) => f.write_str("swap"), } } @@ -196,6 +205,7 @@ where let result = match cli.command { Command::AntiRepeat(args) => anti_repeat::run(args, stdin, stdout, stderr), Command::Bed(args) => bed::run(args, stdin, stdout, stderr), + Command::Compare(args) => compare::run(args, stdin, stdout, stderr), Command::Convert(args) => convert::run(args, stdin, stdout, stderr), Command::Filter(args) => filter::run(args, stdin, stdout, stderr), Command::Liftover(args) => liftover::run(args, stdin, stdout, stderr), @@ -203,6 +213,7 @@ where Command::Score(args) => score::run(args, stdin, stdout, stderr), Command::Split(args) => split::run(args, stdin, stdout, stderr), Command::Sort(args) => sort::run(args, stdin, stdout, stderr), + Command::Stats(args) => stats::run(args, stdin, stdout, stderr), Command::Swap(args) => swap::run(args, stdin, stdout, stderr), }; diff --git a/src/cli/sort.rs b/src/cli/sort.rs index 7f7ae6a..5a8d62a 100644 --- a/src/cli/sort.rs +++ b/src/cli/sort.rs @@ -444,8 +444,12 @@ fn collect_sorted_input( ) -> Result { // Sort reads a single input, so metadata lines are preserved verbatim // (no cross-file deduplication). - let mut accumulator = - SortAccumulator::new(args.sort_by.criterion(), max_in_memory_bytes, temp_dir, false); + let mut accumulator = SortAccumulator::new( + args.sort_by.criterion(), + max_in_memory_bytes, + temp_dir, + false, + ); accumulator.push_stream(reader)?; // Capture counts before `finish` consumes the accumulator. let chains = accumulator.chains_pushed(); @@ -814,7 +818,10 @@ mod tests { String::from_utf8(stdout).unwrap(), "chain 300 chr3 1000 + 0 30 qry1 500 + 0 30 1\n30\n\nchain 100 chr1 1000 + 10 50 qry3 500 + 15 55 2\n40\n\n" ); - assert_eq!(fs::read_to_string(&index.path).unwrap(), "0\t300\n31\t100\n"); + assert_eq!( + fs::read_to_string(&index.path).unwrap(), + "0\t300\n31\t100\n" + ); assert_eq!(stderr, b""); } diff --git a/src/cli/split.rs b/src/cli/split.rs index 41d186f..b95d87b 100644 --- a/src/cli/split.rs +++ b/src/cli/split.rs @@ -801,13 +801,13 @@ impl LoadedInput { let file = File::open(path)?; let mmap = unsafe { memmap2::MmapOptions::new().map(&file)? }; log::debug!("Loaded plain input through memory mapping"); - return Ok(Self { + Ok(Self { bytes: SharedInputBytes::Mmap(Arc::new(mmap)), basename: derive_basename(path), source_path: Some(fs::canonicalize(path)?), source_encoding: SourceEncoding::Plain, _temp: None, - }); + }) } #[cfg(not(feature = "mmap"))] @@ -865,13 +865,13 @@ impl LoadedInput { let file = File::open(&temp.path)?; let mmap = unsafe { memmap2::MmapOptions::new().map(&file)? }; log::debug!("Staged standard input to a temporary mmap-backed file"); - return Ok(Self { + Ok(Self { bytes: SharedInputBytes::Mmap(Arc::new(mmap)), basename: "stdin".to_owned(), source_path: None, source_encoding: SourceEncoding::Plain, _temp: Some(temp), - }); + }) } #[cfg(not(feature = "mmap"))] @@ -1016,7 +1016,7 @@ fn create_temp_path(prefix: &str) -> Result { fn derive_basename(path: &Path) -> String { let file_name = path .file_name() - .unwrap_or_else(|| path.as_os_str()) + .unwrap_or(path.as_os_str()) .to_string_lossy(); let without_gz = file_name.strip_suffix(".gz").unwrap_or(&file_name); let without_chain = without_gz.strip_suffix(".chain").unwrap_or(without_gz); diff --git a/src/cli/stats.rs b/src/cli/stats.rs new file mode 100644 index 0000000..12a8e59 --- /dev/null +++ b/src/cli/stats.rs @@ -0,0 +1,592 @@ +// Copyright (c) 2026 Alejandro Gonzales-Irribarren +// Distributed under the terms of the Apache License, Version 2.0. + +use std::collections::BTreeMap; +use std::io::{BufRead, Write}; +use std::path::PathBuf; + +use chaintools::{ChainError, OwnedChain, Strand, StreamingReader}; +use clap::Args; + +use super::CliError; + +#[derive(Debug, Args)] +pub struct StatsArgs { + #[arg( + short = 'c', + long = "chain", + value_name = "PATH", + help = "Path to the input .chain file. If omitted, read standard input" + )] + chain: Option, + + #[arg(long, help = "Add per-target fragmentation and top-chain tables")] + by_sequence: bool, + + #[arg( + long, + value_name = "N", + default_value_t = 3, + help = "Top scored chains to show per target with --by-sequence" + )] + top: usize, +} + +#[derive(Debug, Clone)] +pub(crate) struct ChainSummary { + pub(crate) id: u64, + pub(crate) score: i64, + pub(crate) aligned_bp: u64, + pub(crate) target_span: u64, + pub(crate) query_span: u64, + pub(crate) blocks: u64, + pub(crate) target_gap_bp: u64, + pub(crate) query_gap_bp: u64, + pub(crate) largest_gap: u64, +} + +#[derive(Debug, Default)] +pub(crate) struct TargetStats { + pub(crate) total_aligned_bp: u64, + pub(crate) top1_share: f64, + pub(crate) top3_share: f64, + pub(crate) chains_50: u64, + pub(crate) chains_90: u64, + pub(crate) chains_95: u64, + pub(crate) chain_n50: u64, + pub(crate) chain_summaries: Vec, +} + +impl TargetStats { + pub(crate) fn chain_count(&self) -> u64 { + self.chain_summaries.len() as u64 + } + + fn finish(&mut self) { + let mut aligned = self + .chain_summaries + .iter() + .map(|chain| chain.aligned_bp) + .collect::>(); + aligned.sort_unstable_by(|a, b| b.cmp(a)); + self.total_aligned_bp = aligned.iter().sum(); + self.top1_share = ratio(aligned.first().copied().unwrap_or(0), self.total_aligned_bp); + self.top3_share = ratio(aligned.iter().take(3).copied().sum(), self.total_aligned_bp); + self.chains_50 = chains_for_fraction(&aligned, self.total_aligned_bp, 50); + self.chains_90 = chains_for_fraction(&aligned, self.total_aligned_bp, 90); + self.chains_95 = chains_for_fraction(&aligned, self.total_aligned_bp, 95); + self.chain_n50 = n50_l50(&aligned).0; + self.chain_summaries.sort_unstable_by(|a, b| { + b.score + .cmp(&a.score) + .then_with(|| b.aligned_bp.cmp(&a.aligned_bp)) + .then_with(|| a.id.cmp(&b.id)) + }); + } +} + +#[derive(Debug, Default)] +pub(crate) struct ChainStats { + pub(crate) chains: u64, + pub(crate) blocks: u64, + pub(crate) aligned_pair_bp: u64, + pub(crate) target_span_bp: u64, + pub(crate) query_span_bp: u64, + pub(crate) target_gap_bp: u64, + pub(crate) query_gap_bp: u64, + pub(crate) target_only_gap_events: u64, + pub(crate) query_only_gap_events: u64, + pub(crate) dual_gap_events: u64, + pub(crate) largest_target_gap: u64, + pub(crate) largest_query_gap: u64, + pub(crate) plus_chains: u64, + pub(crate) minus_chains: u64, + pub(crate) plus_aligned_bp: u64, + pub(crate) minus_aligned_bp: u64, + pub(crate) score_sum: i128, + pub(crate) largest_chain_aligned_bp: u64, + pub(crate) chain_aligned_n50: u64, + pub(crate) chain_aligned_l50: u64, + pub(crate) median_chain_aligned_bp: f64, + pub(crate) median_blocks_per_chain: f64, + pub(crate) median_target_top1_share: f64, + pub(crate) median_target_top3_share: f64, + pub(crate) median_target_chains_90: f64, + pub(crate) max_target_chains_90: u64, + pub(crate) target_sizes: BTreeMap, u32>, + pub(crate) query_sizes: BTreeMap, u32>, + pub(crate) targets: BTreeMap, TargetStats>, + chain_aligned_bps: Vec, + chain_block_counts: Vec, +} + +impl ChainStats { + fn add_chain(&mut self, chain: &OwnedChain) -> Result<(), CliError> { + record_size( + &mut self.target_sizes, + &chain.reference_name, + chain.reference_size, + "target", + )?; + record_size( + &mut self.query_sizes, + &chain.query_name, + chain.query_size, + "query", + )?; + + let aligned_bp = chain.blocks.iter().map(|block| u64::from(block.size)).sum(); + let target_span = u64::from(chain.reference_end - chain.reference_start); + let query_span = u64::from(chain.query_end - chain.query_start); + let mut target_gap_bp = 0u64; + let mut query_gap_bp = 0u64; + let mut largest_gap = 0u64; + + for block in chain + .blocks + .iter() + .take(chain.blocks.len().saturating_sub(1)) + { + let target_gap = u64::from(block.gap_reference); + let query_gap = u64::from(block.gap_query); + target_gap_bp += target_gap; + query_gap_bp += query_gap; + self.largest_target_gap = self.largest_target_gap.max(target_gap); + self.largest_query_gap = self.largest_query_gap.max(query_gap); + largest_gap = largest_gap.max(target_gap).max(query_gap); + match (target_gap > 0, query_gap > 0) { + (true, false) => self.target_only_gap_events += 1, + (false, true) => self.query_only_gap_events += 1, + (true, true) => self.dual_gap_events += 1, + (false, false) => {} + } + } + + self.chains += 1; + self.blocks += chain.blocks.len() as u64; + self.aligned_pair_bp += aligned_bp; + self.target_span_bp += target_span; + self.query_span_bp += query_span; + self.target_gap_bp += target_gap_bp; + self.query_gap_bp += query_gap_bp; + self.score_sum += i128::from(chain.score); + self.largest_chain_aligned_bp = self.largest_chain_aligned_bp.max(aligned_bp); + self.chain_aligned_bps.push(aligned_bp); + self.chain_block_counts.push(chain.blocks.len() as u64); + + match chain.query_strand { + Strand::Plus => { + self.plus_chains += 1; + self.plus_aligned_bp += aligned_bp; + } + Strand::Minus => { + self.minus_chains += 1; + self.minus_aligned_bp += aligned_bp; + } + } + + self.targets + .entry(chain.reference_name.clone()) + .or_default() + .chain_summaries + .push(ChainSummary { + id: chain.id, + score: chain.score, + aligned_bp, + target_span, + query_span, + blocks: chain.blocks.len() as u64, + target_gap_bp, + query_gap_bp, + largest_gap, + }); + Ok(()) + } + + fn finish(&mut self) { + let (n50, l50) = n50_l50(&self.chain_aligned_bps); + self.chain_aligned_n50 = n50; + self.chain_aligned_l50 = l50; + self.median_chain_aligned_bp = median_u64(&self.chain_aligned_bps); + self.median_blocks_per_chain = median_u64(&self.chain_block_counts); + + for target in self.targets.values_mut() { + target.finish(); + } + self.median_target_top1_share = + median_f64(self.targets.values().map(|target| target.top1_share)); + self.median_target_top3_share = + median_f64(self.targets.values().map(|target| target.top3_share)); + self.median_target_chains_90 = median_u64( + &self + .targets + .values() + .map(|target| target.chains_90) + .collect::>(), + ); + self.max_target_chains_90 = self + .targets + .values() + .map(|target| target.chains_90) + .max() + .unwrap_or(0); + } +} + +pub fn run( + args: StatsArgs, + stdin: &mut R, + stdout: &mut W, + _stderr: &mut E, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, + E: Write, +{ + if args.top == 0 { + return Err(CliError::Message( + "--top must be greater than zero".to_owned(), + )); + } + super::ensure_inputs_exist(&[], &[("input chain", args.chain.as_deref())])?; + + let stats = if let Some(path) = &args.chain { + analyze(&mut StreamingReader::from_path(path)?)? + } else { + analyze(&mut StreamingReader::new(stdin))? + }; + write_stats(stdout, &stats, args.by_sequence, args.top)?; + stdout.flush()?; + super::log_summary( + "stats", + &[("chains", stats.chains), ("blocks", stats.blocks)], + ); + Ok(()) +} + +/// Streams one chain input and returns all statistics needed by both commands. +pub(crate) fn analyze(reader: &mut StreamingReader) -> Result { + let mut stats = ChainStats::default(); + while let Some(header) = reader.next_header()? { + let offset = header.offset; + let blocks = reader.read_blocks(offset)?; + let chain = header.into_chain(blocks); + validate_chain(&chain, offset)?; + stats.add_chain(&chain)?; + } + stats.finish(); + Ok(stats) +} + +pub(crate) fn validate_chain(chain: &OwnedChain, offset: usize) -> Result<(), CliError> { + if chain.reference_strand != Strand::Plus { + return Err(format_error( + offset, + "chain target strand must be + for canonical chain coordinates", + )); + } + if chain.reference_start >= chain.reference_end || chain.query_start >= chain.query_end { + return Err(format_error(offset, "chain span is empty or inverted")); + } + if chain.reference_end > chain.reference_size || chain.query_end > chain.query_size { + return Err(format_error(offset, "chain span exceeds its sequence size")); + } + if chain.blocks.iter().any(|block| block.size == 0) { + return Err(format_error( + offset, + "chain contains an empty alignment block", + )); + } + super::validate_block_spans(chain, offset)?; + Ok(()) +} + +fn record_size( + sizes: &mut BTreeMap, u32>, + name: &[u8], + size: u32, + side: &str, +) -> Result<(), CliError> { + if let Some(previous) = sizes.get(name) { + if *previous != size { + return Err(CliError::Message(format!( + "conflicting {side} sizes for {}: {previous} and {size}", + String::from_utf8_lossy(name) + ))); + } + } else { + sizes.insert(name.to_vec(), size); + } + Ok(()) +} + +fn format_error(offset: usize, message: impl Into) -> CliError { + CliError::Chain(ChainError::Format { + offset, + msg: message.into().into(), + }) +} + +fn write_stats( + writer: &mut W, + stats: &ChainStats, + by_sequence: bool, + top: usize, +) -> Result<(), CliError> { + writeln!(writer, "GENERAL")?; + writeln!(writer, "chains\t{}", stats.chains)?; + writeln!(writer, "blocks\t{}", stats.blocks)?; + writeln!(writer, "target_sequences\t{}", stats.target_sizes.len())?; + writeln!(writer, "query_sequences\t{}", stats.query_sizes.len())?; + + writeln!(writer, "\nALIGNMENT")?; + writeln!(writer, "aligned_pair_bp\t{}", stats.aligned_pair_bp)?; + writeln!(writer, "target_span_bp\t{}", stats.target_span_bp)?; + writeln!(writer, "query_span_bp\t{}", stats.query_span_bp)?; + + writeln!(writer, "\nCONTINUITY")?; + writeln!( + writer, + "largest_chain_aligned_bp\t{}", + stats.largest_chain_aligned_bp + )?; + writeln!(writer, "chain_aligned_bp_n50\t{}", stats.chain_aligned_n50)?; + writeln!(writer, "chain_aligned_bp_l50\t{}", stats.chain_aligned_l50)?; + writeln!( + writer, + "median_chain_aligned_bp\t{:.1}", + stats.median_chain_aligned_bp + )?; + writeln!( + writer, + "median_blocks_per_chain\t{:.1}", + stats.median_blocks_per_chain + )?; + writeln!( + writer, + "median_target_top1_share\t{}", + percent(stats.median_target_top1_share) + )?; + writeln!( + writer, + "median_target_top3_share\t{}", + percent(stats.median_target_top3_share) + )?; + writeln!( + writer, + "median_target_chains_to_90pct\t{:.1}", + stats.median_target_chains_90 + )?; + writeln!( + writer, + "max_target_chains_to_90pct\t{}", + stats.max_target_chains_90 + )?; + + writeln!(writer, "\nGAPS")?; + writeln!(writer, "target_gap_bp\t{}", stats.target_gap_bp)?; + writeln!(writer, "query_gap_bp\t{}", stats.query_gap_bp)?; + writeln!( + writer, + "target_only_gap_events\t{}", + stats.target_only_gap_events + )?; + writeln!( + writer, + "query_only_gap_events\t{}", + stats.query_only_gap_events + )?; + writeln!(writer, "dual_gap_events\t{}", stats.dual_gap_events)?; + writeln!(writer, "largest_target_gap\t{}", stats.largest_target_gap)?; + writeln!(writer, "largest_query_gap\t{}", stats.largest_query_gap)?; + + writeln!(writer, "\nSTRAND")?; + writeln!(writer, "plus_chains\t{}", stats.plus_chains)?; + writeln!(writer, "minus_chains\t{}", stats.minus_chains)?; + writeln!(writer, "plus_aligned_bp\t{}", stats.plus_aligned_bp)?; + writeln!(writer, "minus_aligned_bp\t{}", stats.minus_aligned_bp)?; + + writeln!(writer, "\nSCORE")?; + writeln!(writer, "score_sum\t{}", stats.score_sum)?; + writeln!( + writer, + "mean_score_per_chain\t{:.2}", + if stats.chains == 0 { + 0.0 + } else { + stats.score_sum as f64 / stats.chains as f64 + } + )?; + + if by_sequence { + writeln!(writer, "\nBY TARGET")?; + writeln!( + writer, + "target\tchains\taligned_bp\ttop1_share\ttop3_share\tchains_to_50pct\tchains_to_90pct\tchains_to_95pct\tchain_n50" + )?; + for (name, target) in &stats.targets { + writeln!( + writer, + "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}", + String::from_utf8_lossy(name), + target.chain_count(), + target.total_aligned_bp, + percent(target.top1_share), + percent(target.top3_share), + target.chains_50, + target.chains_90, + target.chains_95, + target.chain_n50 + )?; + } + for (name, target) in &stats.targets { + writeln!(writer, "\nTOP CHAINS\t{}", String::from_utf8_lossy(name))?; + write_top_chains(writer, &target.chain_summaries, top)?; + } + } + Ok(()) +} + +pub(crate) fn write_top_chains( + writer: &mut W, + chains: &[ChainSummary], + top: usize, +) -> Result<(), CliError> { + writeln!( + writer, + "rank\tchain_id\tscore\taligned_bp\ttarget_span\tquery_span\tblocks\tdensity\ttarget_gap_bp\tquery_gap_bp\tlargest_gap" + )?; + for (index, chain) in chains.iter().take(top).enumerate() { + writeln!( + writer, + "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}", + index + 1, + chain.id, + chain.score, + chain.aligned_bp, + chain.target_span, + chain.query_span, + chain.blocks, + percent(ratio(chain.aligned_bp, chain.target_span)), + chain.target_gap_bp, + chain.query_gap_bp, + chain.largest_gap + )?; + } + Ok(()) +} + +fn ratio(numerator: u64, denominator: u64) -> f64 { + if denominator == 0 { + 0.0 + } else { + numerator as f64 / denominator as f64 + } +} + +fn percent(value: f64) -> String { + format!("{:.2}%", value * 100.0) +} + +fn n50_l50(values: &[u64]) -> (u64, u64) { + let mut sorted = values.to_vec(); + sorted.sort_unstable_by(|a, b| b.cmp(a)); + let total = sorted.iter().map(|&value| u128::from(value)).sum::(); + let mut cumulative = 0u128; + for (index, value) in sorted.into_iter().enumerate() { + cumulative += u128::from(value); + if cumulative * 2 >= total { + return (value, index as u64 + 1); + } + } + (0, 0) +} + +fn chains_for_fraction(sorted: &[u64], total: u64, percent: u64) -> u64 { + if total == 0 { + return 0; + } + let mut cumulative = 0u128; + for (index, value) in sorted.iter().copied().enumerate() { + cumulative += u128::from(value); + if cumulative * 100 >= u128::from(total) * u128::from(percent) { + return index as u64 + 1; + } + } + sorted.len() as u64 +} + +fn median_u64(values: &[u64]) -> f64 { + let mut sorted = values.to_vec(); + sorted.sort_unstable(); + let middle = sorted.len() / 2; + match sorted.len() { + 0 => 0.0, + length if length % 2 == 1 => sorted[middle] as f64, + _ => (sorted[middle - 1] as f64 + sorted[middle] as f64) / 2.0, + } +} + +fn median_f64(values: impl Iterator) -> f64 { + let mut sorted = values.collect::>(); + sorted.sort_unstable_by(f64::total_cmp); + let middle = sorted.len() / 2; + match sorted.len() { + 0 => 0.0, + length if length % 2 == 1 => sorted[middle], + _ => (sorted[middle - 1] + sorted[middle]) / 2.0, + } +} + +#[cfg(test)] +mod tests { + use super::*; + use std::io::Cursor; + + fn analyze_text(text: &str) -> ChainStats { + analyze(&mut StreamingReader::new(Cursor::new(text.as_bytes()))) + .expect("valid chain fixture") + } + + #[test] + fn collects_alignment_gaps_and_strands() { + let stats = analyze_text( + "chain 100 chr1 1000 + 0 46 qry1 1000 + 0 48 1\n\ + 10 2 0\n10 0 3\n10 4 5\n10\n\n\ + chain 50 chr2 100 + 0 10 qry2 100 - 20 30 2\n10\n\n", + ); + + assert_eq!(stats.chains, 2); + assert_eq!(stats.blocks, 5); + assert_eq!(stats.aligned_pair_bp, 50); + assert_eq!((stats.target_gap_bp, stats.query_gap_bp), (6, 8)); + assert_eq!(stats.target_only_gap_events, 1); + assert_eq!(stats.query_only_gap_events, 1); + assert_eq!(stats.dual_gap_events, 1); + assert_eq!((stats.largest_target_gap, stats.largest_query_gap), (4, 5)); + assert_eq!((stats.plus_chains, stats.minus_chains), (1, 1)); + assert_eq!((stats.plus_aligned_bp, stats.minus_aligned_bp), (40, 10)); + } + + #[test] + fn computes_n50_and_target_concentration_from_aligned_bp() { + let stats = analyze_text( + "chain 4 chrN 1000 + 0 40 q1 1000 + 0 40 1\n40\n\n\ + chain 3 chrN 1000 + 100 130 q1 1000 + 100 130 2\n30\n\n\ + chain 2 chrN 1000 + 200 220 q1 1000 + 200 220 3\n20\n\n\ + chain 1 chrN 1000 + 300 310 q1 1000 + 300 310 4\n10\n\n", + ); + let target = stats.targets.get(b"chrN".as_slice()).expect("target"); + + assert_eq!((stats.chain_aligned_n50, stats.chain_aligned_l50), (30, 2)); + assert_eq!(stats.median_chain_aligned_bp, 25.0); + assert!((target.top1_share - 0.4).abs() < f64::EPSILON); + assert!((target.top3_share - 0.9).abs() < f64::EPSILON); + assert_eq!( + (target.chains_50, target.chains_90, target.chains_95), + (2, 3, 4) + ); + assert_eq!(target.chain_n50, 30); + } +} diff --git a/src/io/index.rs b/src/io/index.rs index 7704e01..4a1f55c 100644 --- a/src/io/index.rs +++ b/src/io/index.rs @@ -121,7 +121,7 @@ impl ChainIndex { { let file = std::fs::File::open(path)?; let mmap = unsafe { MmapOptions::new().map(&file)? }; - return ChainIndex::from_bytes(SharedBytes::from_mmap(mmap)); + ChainIndex::from_bytes(SharedBytes::from_mmap(mmap)) } #[cfg(not(feature = "mmap"))] diff --git a/src/io/reader.rs b/src/io/reader.rs index d053294..2758b60 100644 --- a/src/io/reader.rs +++ b/src/io/reader.rs @@ -238,7 +238,7 @@ impl Reader { { let file = std::fs::File::open(path)?; let mmap = unsafe { MmapOptions::new().map(&file)? }; - return Self::build(SharedBytes::from_mmap(mmap), ParseStrategy::Parallel); + Self::build(SharedBytes::from_mmap(mmap), ParseStrategy::Parallel) } #[cfg(not(feature = "mmap"))] diff --git a/tests/chains.rs b/tests/chains.rs index 5511045..e7b4611 100644 --- a/tests/chains.rs +++ b/tests/chains.rs @@ -12,6 +12,7 @@ fn strand_creation() { } #[test] +#[allow(clippy::clone_on_copy)] fn strand_copy_and_clone() { let original = Strand::Plus; let copied = original;