From 40148c070b21640c2576c55c1f60bde2c5ce3714 Mon Sep 17 00:00:00 2001 From: gwct Date: Thu, 22 Oct 2026 10:15:20 -0400 Subject: [PATCH 1/2] Add per-species scaffolds to stats, progress reporting with --progress-interval, and a fetch BED pre-check --- DEVELOPMENT.md | 118 +++++- README.md | 156 ++++++- mafutils/fetch.py | 369 ++++++++++++++--- mafutils/gc.py | 47 ++- mafutils/index.py | 33 +- mafutils/info.py | 15 +- mafutils/lib/common.py | 181 +++++++- mafutils/lib/options.py | 42 ++ mafutils/stats.py | 410 +++++++++++++++--- tests/example-tabs.maf | 2 +- tests/example-tabs.maf.block.idx | 4 +- tests/example-tabs.maf.scaffold.idx | 4 +- tests/example.maf | 2 +- tests/example.maf.bgz | Bin 338 -> 305 bytes tests/example.maf.bgz.block.idx | 22 +- tests/example.maf.bgz.scaffold.idx | 14 +- tests/example.maf.block.idx | 4 +- tests/example.maf.gz | Bin 285 -> 281 bytes tests/example.maf.gz.block.idx | 4 +- tests/example.maf.gz.scaffold.idx | 4 +- tests/example.maf.scaffold.idx | 4 +- tests/example.maf.scaffold.regenerated.idx | 4 +- tests/expected-maf/08-span-multiple-gaps.maf | 2 +- tests/test_bed_check.py | 163 ++++++++ tests/test_fetch_order.py | 55 +++ tests/test_progress.py | 412 +++++++++++++++++++ tests/test_scaffold_mode.py | 2 +- tests/test_species_scaffolds.py | 273 ++++++++++++ 28 files changed, 2158 insertions(+), 188 deletions(-) create mode 100644 mafutils/lib/options.py create mode 100644 tests/test_bed_check.py create mode 100644 tests/test_fetch_order.py create mode 100644 tests/test_progress.py create mode 100644 tests/test_species_scaffolds.py diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index a0774bd..192dfb1 100644 --- a/DEVELOPMENT.md +++ b/DEVELOPMENT.md @@ -51,6 +51,118 @@ for user-facing installation/usage docs. index time costs 2.06x (0.79s -> 1.62s per 1.5GB), which would undo most of the v0.6.0 speedup. Verified end to end: 1.3s on the 42.7GB file, and the fallback's scanned values equal the header's exactly. +- **Per-species source scaffolds live in `stats`, not `index` or `info`.** + Scaffold names and `srcSize` are only on non-reference `s` lines: parsing + them in `index` measured 2.06x slower indexing, and sampling them in `info` + saturates far too slowly (54 species/scaffold pairs after 1,000 blocks vs + 1,587 after 200,000 and still climbing). `parseBlock` already splits every + `s` line, so `fields[5]` is free there; one `partition(".")` replaces + `speciesFromSrc` and gives an identical species. Keyed by **(species, + scaffold)**: `human.chr4` and `chimp.chr4` are different sequences. + Everything it reports covers scaffolds **present in the MAF**, not the + assembly -- see the README's scope section; keep that caveat (and the + `_in_maf` column names) wherever these numbers surface. +- **Conflicting `srcSize` for one species.scaffold fails the run.** + `reconcileSrcSize` is applied at all three merge levels (lines within a + block, blocks within a chunk, chunks in `mergeStatsInto`) and raises + `SrcSizeConflictError`, which `run_stats` turns into `LOG.error` + exit 1 + without a traceback. Picking one value would be a silent fallback, and + "first seen" would differ across `-p`/`--chunk-size`. A `None` size (did not + parse) never conflicts and is filled in by any line that does parse. Real + data has had zero conflicts (2,647 scaffolds / 500,000 blocks; 1,752 pairs + on the 1.5GB slice), so this is defensive. +- **`tests/example.maf` had an inconsistent `srcSize` from its first commit** + -- `chimp.chr4` was 1000 in one block and 998 in another; nothing read + `srcSize` until the conflict check above. Fixed to 1000 in `example.maf`, + `example-tabs.maf` and `expected-maf/08-span-multiple-gaps.maf`, with the + `.gz`/`.bgz` copies and every index regenerated (only the last block's end + offset moved, +1 byte; `test_scaffold_mode.py`'s hand-written header went + 776 -> 777). Lesson: the script that first "hand-derived" the expected + values assigned each line's size over the previous one, which hid the + conflict. When deriving expected values, collect every distinct value per + key, never overwrite. +- **Progress lines (`COMMON.BlockProgress` / `ByteProgress`)** are logged at + most every `PROGRESS_INTERVAL_SECONDS` (60). `stats` and parallel `gc` call + them once per finished chunk in the main process, against the header's + `blocks=` total. `index` and sequential `gc` see every block, so they check + only every `PROGRESS_CHECK_BLOCKS` (10,000) blocks -- sequential `gc` on + `a` lines only, so `s` lines pay nothing -- and measure bytes via + `diskPosition()`, an `lseek` on the stream's file descriptor. That gives + on-disk (compressed, for gz/bgzip) bytes for every compression type, which + the stream's own offsets (decompressed for gz, virtual for bgzip) do not; + `_HashingRawIO` passes `fileno()` through for it. Measured old -> new, 3 + interleaved runs: `index` on 1.5GB **4.96s -> 4.82s**, sequential `gc` on + 300MB **9.67s -> 9.50s** -- no cost beyond noise -- with byte-identical + outputs. +- **`fetch` progress uses shared memory, not messages** (`newProgressSlots`): + one `multiprocessing.RawArray` int64 per task (batch or scaffold), written + only by that task -- no lock, no contention at any worker count, nothing sent + between processes -- and summed by the main process when its + `wait(timeout=progress.dueIn())` wakes. Needed because block mode has ONE + batch per worker, so the stats-style "update as chunks finish" would say + nothing until the end. Worker cost: one ~45ns store per region (block mode) + or per 8 MB copied (scaffold mode, uncompressed/gzip); bgzip scaffold copies + run per line, so they only count finished scaffolds. Passing a `RawArray` + through the pool initializer was checked under Python 3.14's forkserver + start method. Measured on the full 42.7GB hamster MAF + 979,442-region BED, + `-p 8`, 3 interleaved runs each: committed **101.6s** (100.5-102.8, after + discarding one cold-cache 133.6s run) vs new **102.7s** (101.4-103.7) -- + within run-to-run spread on a shared node whose load varied 11-30. The + computed cost, 979,442 x 45ns = 0.044s of CPU across 8 workers, is far below + what the benchmark can resolve. All 979,442 FASTAs byte-identical. +- **`--progress-interval` is one shared option on each command, not a global + one** (`lib/options.py`, `ProgressIntervalOption`). A global typer option + has to come before the subcommand, so the natural `mafutils stats ... + --progress-interval 300` would fail, and it would be missing from each + command's `--help`. Defining the `Annotated` option once keeps the name, + default, help and validation identical in `index`/`stats`/`gc`/`fetch`. + - **Validation:** 60 (`PROGRESS_INTERVAL_SECONDS`) is the default and the + minimum; `0` is off; 1-59 and negatives are a usage error (exit 2), never + clamped. The minimum is enforced only at the CLI, so `run_*` and + `BlockProgress` accept any interval for tests. + - **Off means nothing runs:** `BlockProgress.enabled` is false, so + `update`/`updateTo` return before reading the clock and `dueIn()` is + `None`, making fetch's `wait()` the plain blocking wait. `index` and + sequential `gc` test a local `report_progress` before even the modulo. + `fetch` allocates no slots, so workers skip the store, and scaffold copies + get no callback. + - **Measured on the 1.5GB slice** (3 interleaved runs each), committed vs + default vs off, outputs byte-identical: `index` 5.10s / 4.85s / 4.85s; + sequential `gc` (300MB) 10.08s / 9.27s / 9.87s. All within noise. +- **`maf_fetch_summary.tsv` is written in batch order.** It used to be written + in the main process as each batch finished, so row order followed worker + timing and changed from run to run (the six benchmark runs above had six + different unsorted md5s, identical sorted). Workers now write their rows to + a per-batch temp file and the parent concatenates them in batch-number order + -- the same approach the loci table already used -- so the file is the same + for any `-p` and in index (genomic) order. Side benefits: the rows are no + longer pickled back through the result pipe, and the worker no longer holds + them in a list. The summary file now appears at the end of the run. +- **Block-mode BED pre-check (`checkBlockRegions`)** runs after the index-order + sort and basename assignment (`-b count` numbers in that order) and before + batching (`chunker` cannot take zero regions: `range()` with step 0). It + reports everything in one error and exits 1 before any output. It replaced + four separate failures: + - `-b id` logged an error per unnamed region, ran anyway, then crashed in a + worker with `KeyError('output_basename')`; + - `parseBed`'s "skip malformed line" passed `file=` to a logger, so the + first malformed line -- or a UCSC `track` line -- crashed with a + `TypeError`; + - inverted or empty intervals were accepted and reported only as "no + overlapping blocks"; + - duplicate output names (the hamster BED has 27,485 repeated 4th-column + values such as `lod=53`) silently overwrote each other's files, which are + opened "w". + On 979,442 regions the check takes 0.2-0.7s and 27-50MB transient. + `--single-output` skips the name checks, since it writes no per-region files. +- **Verified on real data** (1.5GB hamster slice, 258,186 blocks): all 1,752 + species/scaffold rows -- `src_size`, `blocks`, `aligned_bases` -- equal an + independent awk tally exactly, and every pre-existing output (`block.tsv`, + `overall.tsv` apart from the added row, `species.tsv`'s original 13 columns) + is byte-identical to the previous release's. Cost, 3 interleaved runs each, + previous -> new: `-p 1` on a 300MB slice **42.20s -> 42.30s (+0.2%, within + run-to-run spread)**; `-p 8` on the 1.5GB slice **30.44s -> 31.29s + (+2.8%)**, measured while node load varied between 9 and 13. ## Compression / index internals @@ -520,7 +632,11 @@ computed `overall.tsv`/`species.tsv` values against numbers hand-derived directly from `tests/example.maf`'s raw content (not just "doesn't crash"), including a regression guard that the sequential (`-p 1`, no `ProcessPoolExecutor`) and real multi-worker (`-p 2`+) paths produce -identical output. +identical output. `tests/test_species_scaffolds.py` does the same for +`species_scaffolds.tsv` and the two `_in_maf` columns: hand-derived values, +the sum(`aligned_bases`) == `nongaps_total` invariant, invariance to +`-p`/`--chunk-size`, and the conflicting/unparseable `srcSize` paths (run on a +modified copy of `example.maf` in `tmp_path`, never the committed fixture). `tests/test_real_data.py` runs `stats`/`gc`/`fetch` against `tests/real-excerpt.maf` -- 8 real, complete alignment blocks extracted diff --git a/README.md b/README.md index b3e1d6d..871968b 100644 --- a/README.md +++ b/README.md @@ -143,6 +143,21 @@ check for this every run, at two tiers: For a "check once, trust thereafter" workflow instead of passing `--verify-hash` on every invocation, use `mafutils validate` (see below). +## Progress Reporting + +`index`, `stats`, `gc` and `fetch` log a progress line while they run, every +60 seconds by default. Each takes `--progress-interval SECONDS` to change it: + +- `60` or more sets the interval; +- `0` turns progress off entirely: no lines, and none of the bookkeeping + behind them runs; +- anything from 1 to 59, or negative, is rejected with a usage error rather + than adjusted. + +The option belongs to each command, so it goes after the command name: +`mafutils stats input.maf --progress-interval 300`. What the lines count +(blocks, bytes, regions) differs by command; see each command below. + ## Commands ### `mafutils index` @@ -161,6 +176,18 @@ Arguments: | `BLOCK_INDEX` | Output block index path (default: `.block.idx`) | | `SCAFFOLD_INDEX` | Output scaffold index path (default: `.scaffold.idx`) | +Options: + +| Option | Description | +|---|---| +| `--progress-interval` | Seconds between progress lines (default and minimum 60); `0` turns progress off. See Progress Reporting above | + +Every 60 seconds by default it logs how far through the file it is, e.g. `Progress: +12.30 GB / 42.70 GB (28.8%), 123.00 MB/s, 1m 40s elapsed, ~4m 07s left`. +Progress is in bytes because the block count is what indexing produces. For +gzip and bgzip input both numbers are compressed bytes, compared against the +compressed file size. Messages go to the screen only; no log file is written. + ### `mafutils info` Print a quick overview of a MAF without reading it. @@ -211,7 +238,8 @@ index, and collecting them during indexing would mean parsing every non-reference `s` line — measured at 2.06x slower indexing. Sampling the first 1,000 blocks costs milliseconds instead. On a whole-genome MAF those blocks all sit at the start of the first scaffold, so a species absent from that region is -missed. Use `mafutils stats` for an exhaustive per-species breakdown. +missed. Use `mafutils stats` for an exhaustive per-species breakdown, +including each species' source scaffolds and their lengths. **Column meanings.** `ref scaffolds` and `ref bases` are measured on the *reference* sequence. `seq lines` counts every `s` line including the @@ -238,7 +266,7 @@ Options: | Option | Description | |---|---| | `--index`, `-i` | Index file, block-level or scaffold-level matching `--mode`. **Required** — if omitted, looked up at `.block.idx` or `.scaffold.idx` by default; errors if not found there. | -| `--basename`, `-b` | Output basename strategy: `id`, `coords`, or `count` | +| `--basename`, `-b` | Output basename strategy: `id` (the BED's 4th column; every region needs a unique one -- see *BED pre-check* below), `coords`, or `count` | | `--output`, `-o` | Output directory or output filename in single-output mode | | `--fasta`, `-f` | Write FASTA instead of MAF | | `--fasta-header`, `-fh` | FASTA header format: `species-coords-id`, `species-coords`, or `species-only`. The coords-bearing modes always include the source scaffold (`species.scaffold`); `species-only` is deliberately bare so headers match tree tip labels | @@ -252,6 +280,41 @@ Options: | `--verbose` | Emit warning lines from each completed batch | | `--profile` | Log internal timing breakdowns | | `--verify-hash` | Verify the index's stored content hash against the MAF file (see Index Integrity above) | +| `--progress-interval` | Seconds between progress lines (default and minimum 60); `0` turns progress off. See Progress Reporting above | + +**BED pre-check.** In block mode, before doing any work, `fetch` checks the +whole BED and stops with a single error listing every problem it finds: + +- lines it cannot use: fewer than 3 columns, non-integer start or end, a + negative start, or end ≤ start (each reported with its line number); +- with `-b id`, regions that have no name (4th column); +- regions that would write the same output file, where the later one would + silently replace the earlier: names used more than once with `-b id` + (allowed on different scaffolds with `--scaffold-subdirs`), or the same + region listed twice with `-b coords`. `-b count` numbers every region, so it + never collides. + +Blank lines, `#` comments, and UCSC `track`/`browser` header lines are +skipped. Scaffolds that are not in the MAF are not an error; they are reported +as a warning at the end, as before. The check takes under a second on a +979,442-region BED. + +**Progress.** `fetch` logs one progress line every 60 seconds by default +(`--progress-interval`), combining all workers: + +- **Block mode** counts regions, e.g. `Progress: 412,000 / 979,442 regions + (42.1%), 1,850 regions/s, 3m 42s elapsed, ~5m 06s left`. Workers finish + their shares at different times (the "Completed batch N" lines show which), + so the combined rate drops near the end and the last estimates run low. +- **Scaffold mode** counts bytes copied against the total size of the requested + scaffolds, and logs `Finished (k of M scaffolds)` as each one + completes. For bgzip input it counts finished scaffolds instead of bytes: + its copy loop runs once per line and is kept free of progress checks. + +Each worker records its count in its own slot of shared memory, and the main +process reads them only when a line is due. No messages pass between processes +and no worker ever waits on another, so this adds no measurable time at any +scale. ### Reference-side columns in `maf_fetch_summary.tsv` @@ -463,9 +526,97 @@ Options: | `--expected-species-file` | File with one species name per line | | `--verify-hash` | Verify the index's stored content hash against the MAF file (see Index Integrity above) | | `--log-level` | Logging level: `DEBUG`, `INFO`, `WARNING`, `ERROR` | +| `--progress-interval` | Seconds between progress lines (default and minimum 60); `0` turns progress off. See Progress Reporting above | | `--html-dashboard` | Write an HTML summary dashboard | | `--dashboard-top-species` | Number of species shown in dashboard bar plots | | `--dashboard-max-block-points` | Maximum block rows sampled for dashboard plots | +| `--dashboard-top-scaffolds` | Largest scaffolds listed per species in the dashboard (default 10). The TSV is never capped. | + +Outputs: + +| File | Contents | +|---|---| +| `.overall.tsv` | One metric per line | +| `.species.tsv` | One row per species | +| `.species_scaffolds.tsv` | One row per (species, source scaffold) present in the MAF — see below | +| `.block.tsv` | One row per block (skip with `--no-block-table`) | +| `.dashboard.html` | Self-contained summary dashboard (with `--html-dashboard`) | + +**Progress.** While it runs, `stats` logs a progress line every 60 seconds by +default (`--progress-interval`), e.g. `Progress: 1,215,000 / 8,094,203 blocks (15.0%), 2,850 blocks/s, 7m 06s +elapsed, ~40m 13s left`. The total comes from the index header, so an index +built before mafutils 0.7 gives counts and rate only (rebuild it to get the +percentage). Block sizes vary across a genome, so the time estimate is rough. +If the number of blocks processed differs from the header's count, `stats` +warns that the index rows and header disagree and the results may not cover +the whole MAF. `gc` reports progress the same way; `index` measures it in +bytes instead (see `mafutils index`). + +### Per-species source scaffolds + +`stats` reads every `s` line of every block, so it also records which source +scaffolds each species contributes and how long they are. Nothing else reports +this for non-reference species: the index never parses non-reference lines, and +`info` only samples species names. + +**Scope: scaffolds present in the MAF, not the whole assembly.** A scaffold +appears only if at least one of its bases aligned to the reference. Unaligned +scaffolds — unplaced contigs, sequence outside the aligned regions — never +appear in the MAF, so `stats` cannot see them. Each row's length is exact (it is +the MAF's `srcSize`, which the aligner takes from the assembly), but the *set* +of scaffolds is limited to those in the MAF. The per-species totals therefore +err in opposite directions depending on the comparison: + +| Compared with | Direction | Why | +|---|---|---| +| the true assembly size | **under**counts | unaligned scaffolds contribute nothing | +| aligned bases (`nongaps_total`) | **over**counts | a barely-aligned 200 Mb scaffold contributes all 200 Mb | + +The gap to the true assembly cannot be measured from the MAF alone — it needs +the assembly's own scaffold list. It is small for chromosome-level assemblies +and can be large for fragmented ones. Do not read `scaffold_bases_in_maf` as an +assembly-size estimate. + +`.species_scaffolds.tsv`, sorted by species, then largest scaffold +first: + +| Column | Meaning | +|---|---| +| `species` | Species (the part of the MAF `src` before the first `.`) | +| `scaffold` | Source scaffold (the part after it); `.` if the `src` has no `.` | +| `src_size` | Full length of that scaffold in that species' assembly (MAF `srcSize`); `.` if it did not parse as an integer | +| `blocks` | Blocks in which this scaffold appears | +| `aligned_bases` | Non-gap bases from this scaffold across all blocks | +| `aligned_fraction` | `aligned_bases / src_size`. **Not capped at 1.0**: a region aligned into several blocks (duplicated or paralogous alignments) can push it past 1, and that is worth seeing | + +Rows are keyed by (species, scaffold), never by scaffold name alone — +`human.chr4` and `chimp.chr4` are different sequences that happen to share a +name. + +Two columns are appended to the **end** of `.species.tsv`, so no +existing column moves: + +| Column | Meaning | +|---|---| +| `scaffolds_in_maf` | Distinct source scaffolds of this species present in the MAF | +| `scaffold_bases_in_maf` | Sum of those scaffolds' full lengths | + +Per species, the `aligned_bases` of its scaffolds sum exactly to its +`nongaps_total`, and their `blocks` sum to at least its `blocks_present` +(more when a block holds several copies of the species). + +The dashboard adds both columns to its species table and a collapsed section +per species listing the largest scaffolds. Each section's summary line reads +"showing N of M scaffolds", so a truncated list is visible before expanding it. + +**`srcSize` integrity.** If one species.scaffold reports two different +`srcSize` values, the MAF is internally inconsistent and `stats` exits with an +error naming the scaffold and both values — picking one would hide the problem. +An `s` line whose `srcSize` is not an integer still counts toward every other +statistic. Its scaffold takes its length from any other line of that scaffold +that does parse, and is written with `src_size` `.` only if none does. The +number of such lines is logged as a warning and written to `overall.tsv` as +`scaffold_size_parse_errors`. ### `mafutils gc` @@ -497,6 +648,7 @@ Options: | `--chunk-size` | Blocks per worker task in parallel mode | | `--verify-hash` | Verify the index's stored content hash against the MAF file (see Index Integrity above; only checked when an index is actually used) | | `--log-level` | Logging level: `DEBUG`, `INFO`, `WARNING`, `ERROR` | +| `--progress-interval` | Seconds between progress lines (default and minimum 60); `0` turns progress off. See Progress Reporting above | **Sequential vs. parallel processing:** GC counting is cheap per base, so the bottleneck is I/O/decompression rather than CPU. There are two processing diff --git a/mafutils/fetch.py b/mafutils/fetch.py index 5f7aaac..086870e 100644 --- a/mafutils/fetch.py +++ b/mafutils/fetch.py @@ -70,6 +70,9 @@ import statistics import tempfile import traceback +import shutil +import ctypes +import multiprocessing from enum import Enum from types import SimpleNamespace from typing import Annotated, Optional @@ -80,6 +83,7 @@ from mafutils.lib import common as COMMON from mafutils.lib import loginit as LOGINIT +from mafutils.lib import options as OPTIONS ############################################################################# @@ -105,6 +109,9 @@ WORKER_PROFILE = None WORKER_MAF_FP = None WORKER_BLOCK_CACHE = None +# One shared int64 per task, written only by that task -- see newProgressSlots. +WORKER_PROGRESS_SLOTS = None +WORKER_SUMMARY_TMP_DIR = None ############################################################################# @@ -185,6 +192,8 @@ def initBatchWorker( index_file=None, scaffold_runs=None, loci_tmp_dir=None, + progress_slots=None, + summary_tmp_dir=None, ): global WORKER_HEADER global WORKER_MAF_FILE @@ -204,6 +213,8 @@ def initBatchWorker( global WORKER_PROFILE global WORKER_MAF_FP global WORKER_BLOCK_CACHE + global WORKER_PROGRESS_SLOTS + global WORKER_SUMMARY_TMP_DIR WORKER_HEADER = header WORKER_MAF_FILE = maf_file @@ -224,6 +235,8 @@ def initBatchWorker( WORKER_FASTA_DEDUPE = fasta_dedupe WORKER_VERBOSE = verbose WORKER_PROFILE = profile + WORKER_PROGRESS_SLOTS = progress_slots + WORKER_SUMMARY_TMP_DIR = summary_tmp_dir # For gzip input, prefetched_cache already holds every block needed across # all regions, decoded in strictly-ascending file order (see @@ -237,6 +250,24 @@ def initBatchWorker( ############################################################################# +def newProgressSlots(n_tasks): + """ + Shared memory for worker progress: one int64 per task, which only that + task ever writes (a plain store, ~45ns) and the main process sums when a + progress line is due. One writer per slot means no lock and no contention + at any worker count, and nothing is ever sent between processes -- so the + cost per region is the store alone, independent of MAF size. A Manager + value or a queue would instead make every update a round trip to another + process. + """ + return multiprocessing.RawArray(ctypes.c_int64, n_tasks) + + +def initScaffoldWorker(progress_slots): + global WORKER_PROGRESS_SLOTS + WORKER_PROGRESS_SLOTS = progress_slots + + def closeBatchWorker(): global WORKER_MAF_FP if WORKER_MAF_FP is not None: @@ -454,30 +485,113 @@ def parseBed(bed_file, LOG, mode="block"): - basename (from column 4, if present; otherwise None) In "scaffold" mode: accepts 1+ columns; only uses first column (scaffold) + + Returns (regions, problems). `problems` lists (line_number, description) + for block-mode lines that cannot be used -- malformed, or not a valid + interval -- for the caller's BED pre-check to report all at once. They + used to be "skipped with a warning", but that call passed file= to a + logger, so the first bad line crashed fetch with a TypeError instead; + and an inverted or empty interval was accepted and only ever showed up + as "no overlapping blocks". Blank lines, # comments and UCSC + track/browser header lines are skipped. """ regions = [] + problems = [] with open(bed_file, "r") as f: - for line in f: + for line_number, line in enumerate(f, start=1): line = line.strip() if not line or line.startswith("#"): continue fields = line.split() + if fields[0] in ("track", "browser"): + continue if mode == "scaffold": scaffold = fields[0] region = {"scaffold": scaffold, "start": 0, "end": 0, "id": fields[3] if len(fields) >= 4 else None} regions.append(region) - else: - # block mode: require at least 3 columns (scaffold/start/end) - try: - scaffold = fields[0] - start = int(fields[1]) - end = int(fields[2]) - region = {"scaffold": scaffold, "start": start, "end": end, "id": fields[3] if len(fields) >= 4 else None} - regions.append(region) - except (IndexError, ValueError): - LOG.warning(f"Skipping malformed BED line: {line}", file=sys.stderr) - continue - return regions + continue + if len(fields) < 3: + problems.append((line_number, "fewer than 3 columns")) + continue + try: + start = int(fields[1]) + end = int(fields[2]) + except ValueError: + problems.append((line_number, "start/end not integers")) + continue + if start < 0 or end <= start: + problems.append((line_number, f"invalid interval {fields[0]}:{start}-{end}")) + continue + regions.append({"scaffold": fields[0], "start": start, "end": end, "id": fields[3] if len(fields) >= 4 else None}) + return regions, problems + + +def _count(n, noun): + return f"{n:,} {noun}" if n == 1 else f"{n:,} {noun}s" + + +def checkBlockRegions(regions, bed_problems, basename_mode, single_output, scaffold_subdirs): + """ + The block-mode BED pre-check: every problem that would make a run crash + part-way or silently lose output, found before any work and reported + together. Returns one short message per kind of problem (empty if none), + each with a single example. Expects output_basename already assigned. + + Collisions are what matter for output: two regions whose output files have + the same path means the later silently replaces the earlier (files are + opened "w"). With --scaffold-subdirs the scaffold is part of the path, so + the same name on different scaffolds is fine. Memory is one set entry per + region, freed on return. + """ + messages = [] + if bed_problems: + line_number, why = bed_problems[0] + messages.append(f"{_count(len(bed_problems), 'unusable line')}, e.g. line {line_number}: {why}") + + if not regions: + messages.append("no usable regions") + return messages + + if single_output: + return messages # one combined output: no per-region names + + if basename_mode == "id": + unnamed = [r for r in regions if not r.get("id")] + if len(unnamed) == len(regions): + messages.append("no region has a name (4th column), which -b id needs") + elif unnamed: + r = unnamed[0] + messages.append( + f"{len(unnamed):,} of {_count(len(regions), 'region')} have no name (4th column), which -b id " + f"needs, e.g. {r['scaffold']}:{r['start']}-{r['end']}" + ) + + seen = defaultdict(set) # scaffold (or "" without subdirs) -> names + repeats = defaultdict(int) # (scaffold or "", name) -> extra occurrences + for r in regions: + name = r.get("output_basename") + if name is None: # unnamed under -b id: reported above + continue + scope = r["scaffold"] if scaffold_subdirs else "" + if name in seen[scope]: + repeats[(scope, name)] += 1 + else: + seen[scope].add(name) + if repeats: + (_scope, name), extra = next(iter(repeats.items())) + verb = "is" if len(repeats) == 1 else "are" + if basename_mode == "id": + where = " on the same scaffold" if scaffold_subdirs else "" + messages.append( + f"{_count(len(repeats), 'name')} {verb} shared by several regions{where}, whose files would " + f"overwrite each other, e.g. {name} ({extra + 1:,} regions)" + ) + else: # coords: scaffold-start-end, so only identical regions collide + messages.append( + f"{_count(len(repeats), 'region')} {verb} listed more than once, so their files would overwrite " + f"each other, e.g. {'{}:{}-{}'.format(*name.rsplit('-', 2))} ({extra + 1:,} times)" + ) + return messages ############################################################################# @@ -1552,7 +1666,20 @@ def fetchByBatch( BATCHLOG.error(f"Batch {batch_num}: Worker MAF handle is not initialized.") raise RuntimeError("Worker MAF handle is not initialized.") - batch_results = [None] * len(batch) + # Summary rows go to a per-batch file that the parent concatenates in batch + # order (same as the loci table below), so maf_fetch_summary.tsv no longer + # depends on which worker finishes first. Only --single-output, which has + # no summary rows to order, still returns them in memory. + batch_results = [None] * len(batch) if WORKER_SUMMARY_TMP_DIR is None else None + summary_path = None + summary_fp = None + if WORKER_SUMMARY_TMP_DIR is not None: + summary_path = os.path.join(WORKER_SUMMARY_TMP_DIR, f"summary.task{batch_num}.tsv") + summary_fp = open(summary_path, "w", encoding="utf-8") + # This batch's progress slot (see newProgressSlots), as locals so the + # per-region store is one plain write; None with --progress-interval 0. + progress_slots = WORKER_PROGRESS_SLOTS + progress_slot = batch_num - 1 region_num = 0 no_output_regions = 0 written_regions = 0 @@ -1630,7 +1757,12 @@ def fetchByBatch( # conservative) offsets are kept -- correct, just more rescanning. resume_by_scaffold[scaffold] = scan_state["resume"] - batch_results[result_idx] = region_result + if summary_fp is not None: + summary_fp.write(region_result + "\n") + else: + batch_results[result_idx] = region_result + if progress_slots is not None: # None: --progress-interval 0 + progress_slots[progress_slot] = region_num if not WORKER_SINGLE_OUTPUT: # Field 4 is ref.n.overlapping.blocks (see summary_headers). This # used to read field 1, which is the region's START coordinate @@ -1651,13 +1783,16 @@ def fetchByBatch( index_fp.close() if loci_fp is not None: loci_fp.close() + if summary_fp is not None: + summary_fp.close() if profile_state is not None: profile_state["time_batch_total"] = time.perf_counter() - batch_timer_start return { "results": batch_results, - "processed_regions": len(batch_results), + "summary_path": summary_path, + "processed_regions": len(batch), "zero_overlap_regions": no_output_regions, "no_output_regions": no_output_regions, "batch_num": batch_num, @@ -1671,21 +1806,32 @@ def fetchByBatch( ############################################################################# -def fetchByScaffold(scaffold, start_end, maf_file, maf_header, maf_compression, out_dir, LOG, scaffold_subdirs=False): +def fetchByScaffold(scaffold, start_end, maf_file, maf_header, maf_compression, out_dir, LOG, scaffold_subdirs=False, + progress_slot=None): """ Worker function to extract one scaffold from the MAF file. start_end is a tuple (start_byte, end_byte). + + progress_slot: this task's index into WORKER_PROGRESS_SLOTS, which gets + the bytes copied so far after every 8 MB chunk (see newProgressSlots). """ start, end = start_end output_path = buildOutputPath(out_dir, scaffold, scaffold, ".maf", scaffold_subdirs) + on_progress = None + if progress_slot is not None: + slots = WORKER_PROGRESS_SLOTS + + def on_progress(written): + slots[progress_slot] = written + try: with COMMON.openMaf(maf_file, maf_compression, "rb") as mfp, open(output_path, "wb") as outfp: LOG.info(f">>> Scaffold {scaffold}"); outfp.write(maf_header.encode("utf-8")) # Streamed in bounded chunks: a scaffold can be hundreds of GB, # far too large to materialize as one bytes object. - written = COMMON.copyMafRangeToStream(mfp, maf_compression, start, end, outfp) + written = COMMON.copyMafRangeToStream(mfp, maf_compression, start, end, outfp, on_progress=on_progress) except Exception as e: # repr(), plus the traceback, because some of the exceptions that # actually occur here have an empty str() -- MemoryError being the @@ -1697,29 +1843,40 @@ def fetchByScaffold(scaffold, start_end, maf_file, maf_header, maf_compression, return {"scaffold": scaffold, "output_path": output_path, "ok": True, "bytes_written": written} -def fetchScaffoldsSequential(ordered_scaffolds, index, maf_file, maf_header, maf_compression, out_dir, LOG, scaffold_subdirs=False): +def fetchScaffoldsSequential(ordered_scaffolds, index, maf_file, maf_header, maf_compression, out_dir, LOG, scaffold_subdirs=False, + progress=None): """ Extracts multiple scaffolds using a single open MAF handle, in the given (ascending file-offset) order. Used for gzip input: reusing one handle and never seeking backward avoids the redundant-decompression cost that reopening the file per scaffold (as fetchByScaffold does) would incur. + + Runs in the main process, so `progress` (a ByteProgress over the summed + scaffold spans) is fed directly, once per 8 MB chunk. """ results = [] + copied_before = 0 # bytes from scaffolds already finished with COMMON.openMaf(maf_file, maf_compression, "rb") as mfp: - for scaffold in ordered_scaffolds: + for n_done, scaffold in enumerate(ordered_scaffolds, start=1): start, end = index[scaffold] output_path = buildOutputPath(out_dir, scaffold, scaffold, ".maf", scaffold_subdirs) + on_progress = None + if progress is not None: + def on_progress(written, base=copied_before): + progress.updateTo(base + written) try: LOG.info(f">>> Scaffold {scaffold}"); with open(output_path, "wb") as outfp: outfp.write(maf_header.encode("utf-8")) # Streamed in bounded chunks -- see fetchByScaffold. - written = COMMON.copyMafRangeToStream(mfp, maf_compression, start, end, outfp) + written = COMMON.copyMafRangeToStream(mfp, maf_compression, start, end, outfp, on_progress=on_progress) except Exception as e: LOG.error(f"{scaffold}: FAILED to extract: {e!r}") LOG.debug("Traceback for %s:\n%s", scaffold, traceback.format_exc()) results.append({"scaffold": scaffold, "output_path": output_path, "ok": False, "error": repr(e)}) continue + copied_before += written + LOG.info(f"Finished {scaffold} ({n_done} of {len(ordered_scaffolds)} scaffolds)") results.append({"scaffold": scaffold, "output_path": output_path, "ok": True, "bytes_written": written}) return results @@ -1742,11 +1899,12 @@ def run_fetch(args, cmdline="mafutils fetch"): LOG.info(f"mafutils fetch called as: {cmdline}"); LOG.info("-" * 40); + arg_width = max(len(arg) for arg in vars(args)) for arg, value in vars(args).items(): valstr = str(value) if len(valstr) > 80: valstr = valstr[:77] + "..." - LOG.info(f"{arg:20} : {valstr}") + LOG.info(f"{arg:{arg_width}} : {valstr}") LOG.info("-" * 40 ) # Log the command line arguments @@ -1848,7 +2006,7 @@ def run_fetch(args, cmdline="mafutils fetch"): ) LOG.info(f"Parsing BED file...... {args.bed_file}"); - regions = parseBed(args.bed_file, LOG, args.mode); + regions, bed_problems = parseBed(args.bed_file, LOG, args.mode); expected_species = parseExpectedSpecies(args) if expected_species: @@ -1895,6 +2053,26 @@ def run_fetch(args, cmdline="mafutils fetch"): # Process scaffolds in ascending file-offset order; for gzip this is # what lets a single reused handle avoid ever seeking backward. ordered_scaffolds = sorted(scaffold_set, key=lambda s: index[s][0]) + n_scaffolds = len(ordered_scaffolds) + + # Progress in bytes copied, against the summed scaffold spans -- exact + # for uncompressed (file bytes) and gzip (decompressed offsets). bgzip + # offsets are virtual and its copy loop runs per LINE, which is kept + # free of checks, so bgzip counts finished scaffolds instead. + # With --progress-interval 0 there are no slots and no copy callback. + progress_on = args.progress_interval > 0 + byte_progress = progress_on and maf_compression != "bgzip" + LOG.info(COMMON.progressStartupLine(args.progress_interval)) + if byte_progress: + progress = COMMON.ByteProgress(LOG, sum(index[s][1] - index[s][0] for s in ordered_scaffolds), + interval=args.progress_interval) + else: + if progress_on: + LOG.info( + "bgzip input: progress counts finished scaffolds, not bytes (its copy loop runs " + "once per line and is kept free of progress checks)." + ) + progress = COMMON.BlockProgress(LOG, n_scaffolds, interval=args.progress_interval, unit="scaffolds") scaffold_results = [] if maf_compression == "gz": @@ -1904,14 +2082,21 @@ def run_fetch(args, cmdline="mafutils fetch"): "gzip-compressed files; --processes will be ignored. (Use a bgzip-compressed MAF " "for real parallel speedup on compressed input.)" ) - scaffold_results = fetchScaffoldsSequential(ordered_scaffolds, index, args.maf_file, maf_header, maf_compression, args.output, LOG, scaffold_subdirs=args.scaffold_subdirs) + scaffold_results = fetchScaffoldsSequential(ordered_scaffolds, index, args.maf_file, maf_header, maf_compression, args.output, LOG, scaffold_subdirs=args.scaffold_subdirs, progress=progress if byte_progress else None) else: - # Extract from MAF using the region index for matching scaffolds - with ProcessPoolExecutor(max_workers=args.processes) as executor: - futures = [] - for scaffold in ordered_scaffolds: + # Extract from MAF using the region index for matching scaffolds. + # Scaffold i writes its bytes-copied count to slot i (see + # newProgressSlots); results are stored by slot so the reporting + # below sees them in the original order, whatever order they finish. + progress_slots = newProgressSlots(n_scaffolds) if byte_progress else None + results_by_slot = [None] * n_scaffolds + n_finished = 0 + with ProcessPoolExecutor(max_workers=args.processes, initializer=initScaffoldWorker, + initargs=(progress_slots,)) as executor: + futures = {} + for slot, scaffold in enumerate(ordered_scaffolds): start, end = index[scaffold] - futures.append(executor.submit( + futures[executor.submit( fetchByScaffold, scaffold, (start, end), @@ -1921,9 +2106,20 @@ def run_fetch(args, cmdline="mafutils fetch"): args.output, LOG, args.scaffold_subdirs, - )) - for future in futures: - scaffold_results.append(future.result()) + slot if byte_progress else None, + )] = slot + while futures: + # Times out exactly when a progress line is due. + done, _pending = wait(list(futures.keys()), timeout=progress.dueIn(), return_when=FIRST_COMPLETED) + for future in done: + slot = futures.pop(future) + result = future.result() + results_by_slot[slot] = result + n_finished += 1 + status = "Finished" if result["ok"] else "FAILED" + LOG.info(f"{status} {result['scaffold']} ({n_finished} of {n_scaffolds} scaffolds)") + progress.updateTo(sum(progress_slots) if byte_progress else n_finished) + scaffold_results = results_by_slot # Only report "Wrote" for extractions that actually produced data -- # never for one that raised or came back empty. @@ -1998,6 +2194,37 @@ def run_fetch(args, cmdline="mafutils fetch"): unknown_rank = len(scaffold_first_run) regions.sort(key=lambda r: (scaffold_first_run.get(r["scaffold"], unknown_rank), r["start"], r["end"])) + # Pre-assign output basenames for regions -- after the sort, because + # -b count numbers regions per scaffold in this (index) order. + region_counter = {} + for region in regions: + scaffold = region["scaffold"] + start = region["start"] + end = region["end"] + basename_mode = args.basename + + if basename_mode == "id": + # None when the BED has no name; the pre-check below reports that. + region["output_basename"] = region["id"] + elif basename_mode == "count": + count = region_counter.get(scaffold, 0) + 1 + region_counter[scaffold] = count + region["output_basename"] = f"{scaffold}.{count}" + else: # basename_mode == "coords": + region["output_basename"] = f"{scaffold}-{start}-{end}" + + # BED pre-check: everything that would crash the run part-way or silently + # lose output, reported in ONE error before any work (and before batching, + # which cannot handle zero regions). See checkBlockRegions. + bed_check = checkBlockRegions(regions, bed_problems, args.basename, args.single_output, args.scaffold_subdirs) + if len(bed_check) == 1: + LOG.error(f"BED check failed for {args.bed_file}: {bed_check[0]}. Nothing was fetched.") + sys.exit(1) + if bed_check: + LOG.error(f"BED check failed for {args.bed_file}; nothing was fetched:\n - " + "\n - ".join(bed_check)) + sys.exit(1) + LOG.info(f"BED check passed: {_count(len(regions), 'region')}.") + num_regions = len(regions); # One batch per worker: each streams its own span of the index. Splitting # by region count (not by scaffold) keeps parallelism independent of how @@ -2022,27 +2249,6 @@ def run_fetch(args, cmdline="mafutils fetch"): LOG.info(f"Using single output file: {output_file}") ## Single output mode being developed ! - # Pre-assign output basenames for regions. - # For regions missing a 4th column (basename), assign a unique counter per scaffold. - region_counter = {} - for region in regions: - scaffold = region["scaffold"] - start = region["start"] - end = region["end"] - basename_mode = args.basename - - if basename_mode == "id": - if region.get("id"): - region["output_basename"] = region["id"] - else: - LOG.error(f"Basename mode set to 'id', but '{region}' lacks an ID (4th column)."); - elif basename_mode == "count": - count = region_counter.get(scaffold, 0) + 1 - region_counter[scaffold] = count - region["output_basename"] = f"{scaffold}.{count}" - else: # basename_mode == "coords": - region["output_basename"] = f"{scaffold}-{start}-{end}" - if args.scaffold_subdirs: # Pre-create every needed scaffold subdirectory upfront, in the main # process, mirroring the single os.makedirs(args.output, ...) above -- @@ -2070,6 +2276,24 @@ def run_fetch(args, cmdline="mafutils fetch"): loci_tmp = tempfile.TemporaryDirectory(prefix="maf_fetch_loci_", dir=args.output) if args.loci_table else None loci_tmp_dir = loci_tmp.name if loci_tmp is not None else None loci_paths = {} + # Workers write summary rows per batch; concatenated in batch order below. + summary_tmp = None if args.single_output else tempfile.TemporaryDirectory(prefix="maf_fetch_summary_", dir=args.output) + summary_tmp_dir = summary_tmp.name if summary_tmp is not None else None + summary_paths = {} + # ref.* prefix marks the columns measured on the REFERENCE, so they are + # not mistaken for per-species facts: the reference is contiguous by + # construction while a species can jump kilobases at the same boundary + # (that asymmetry is what the n.*/split.* columns exist to report). + # ref.block.bases was "block.lengths", whose plural implied a list when + # it is a single sum; it falls below (end - start) exactly when the + # region has reference bases no block covers. + summary_headers = ["scaffold", "start", "end", "basename", "ref.n.overlapping.blocks", "ref.block.bases", "n.sequences"] + ELEMENT_LOCI_HEADERS + # Batch b writes its regions-done count to slot b-1 (see newProgressSlots); + # the wait() below wakes when a line is due and sums them. + # With --progress-interval 0 there are no slots, so workers skip the store. + progress = COMMON.BlockProgress(LOG, num_regions, interval=args.progress_interval, unit="regions") + LOG.info(COMMON.progressStartupLine(args.progress_interval)) + progress_slots = newProgressSlots(len(batches)) if progress.enabled else None with ProcessPoolExecutor( max_workers=effective_processes, initializer=initBatchWorker, @@ -2091,18 +2315,10 @@ def run_fetch(args, cmdline="mafutils fetch"): index_file, scaffold_runs, loci_tmp_dir, + progress_slots, + summary_tmp_dir, ), - ) as executor, open(info_outfile, "w") as info_out: - # ref.* prefix marks the columns measured on the REFERENCE, so they are - # not mistaken for per-species facts: the reference is contiguous by - # construction while a species can jump kilobases at the same boundary - # (that asymmetry is what the n.*/split.* columns exist to report). - # ref.block.bases was "block.lengths", whose plural implied a list when - # it is a single sum; it falls below (end - start) exactly when the - # region has reference bases no block covers. - summary_headers = ["scaffold", "start", "end", "basename", "ref.n.overlapping.blocks", "ref.block.bases", "n.sequences"] + ELEMENT_LOCI_HEADERS - info_out.write("\t".join(summary_headers) + "\n") - + ) as executor: futures = {} for batch_num, batch in enumerate(batches, start=1): future = executor.submit( @@ -2113,7 +2329,11 @@ def run_fetch(args, cmdline="mafutils fetch"): futures[future] = batch_num while futures: - done, _pending = wait(list(futures.keys()), return_when=FIRST_COMPLETED) + # Times out exactly when a progress line is due, so the main + # process sleeps between lines instead of polling. + done, _pending = wait(list(futures.keys()), timeout=progress.dueIn(), return_when=FIRST_COMPLETED) + if progress_slots is not None: + progress.updateTo(sum(progress_slots)) for future in done: batch_num = futures.pop(future) result = future.result() @@ -2126,8 +2346,7 @@ def run_fetch(args, cmdline="mafutils fetch"): loci_paths[result["batch_num"]] = result["loci_path"] if not args.single_output: - for r in result["results"]: - info_out.write(r + "\n") + summary_paths[result["batch_num"]] = result["summary_path"] total_regions_reported += result["processed_regions"] zero_overlap_regions += result["zero_overlap_regions"] missing_scaffold_regions += result["warning_counts"].get("missing-scaffold", 0) @@ -2184,6 +2403,18 @@ def run_fetch(args, cmdline="mafutils fetch"): LOG.error(fail_fast_message) sys.exit(2) + # In BATCH-NUMBER order, not completion order: batches partition the + # index-ordered regions, so the table is deterministic and in genomic order. + # It used to be written as batches finished, so its row order changed from + # run to run with worker timing. Same approach as the loci table below. + with open(info_outfile, "w", encoding="utf-8") as info_out: + info_out.write("\t".join(summary_headers) + "\n") + for batch_num in sorted(summary_paths): + with open(summary_paths[batch_num], "r", encoding="utf-8") as fp: + shutil.copyfileobj(fp, info_out) + if summary_tmp is not None: + summary_tmp.cleanup() + if args.profile and profile_totals is not None and profiled_batches > 0: LOG.info("Profile summary across %d batch(es):", profiled_batches) LOG.info( @@ -2299,7 +2530,7 @@ def fetch_command( help="Index file (block-level or scaffold-level, matching --mode). Required -- if omitted, looked up at .block.idx or .scaffold.idx by default; errors if not found there.", ), ] = None, - basename: Annotated[BasenameMode, typer.Option("--basename", "-b", help="Basename strategy: 'id' (BED 4th col if present), 'coords' (scaffold-start-end), or 'count' (numbered per scaffold)")] = BasenameMode.coords, + basename: Annotated[BasenameMode, typer.Option("--basename", "-b", help="Basename strategy: 'id' (BED 4th column; every region must have one), 'coords' (scaffold-start-end), or 'count' (numbered per scaffold)")] = BasenameMode.coords, output: Annotated[str, typer.Option("--output", "-o", help="Output directory (default: current directory); if --single-output is used, this should instead be a filename and will default to maf-fetch.maf in the current directory.")] = ".", fasta: Annotated[bool, typer.Option("--fasta", "-f", help="Output region(s) in FASTA format instead of MAF.")] = False, fasta_header: Annotated[FastaHeaderMode, typer.Option("--fasta-header", "-fh", help="Controls the FASTA header format: 'species-coords-id' (default, includes species, coordinates, and region ID if present),'species-coords' (species and coordinates only), or 'species-only' (just the species name).")] = FastaHeaderMode.species_coords, @@ -2314,6 +2545,7 @@ def fetch_command( verbose: Annotated[bool, typer.Option("--verbose", help="Print every warning line after each batch completes.")] = False, profile: Annotated[bool, typer.Option("--profile", help="Log aggregate timing breakdowns for internal fetch steps.")] = False, verify_hash: Annotated[bool, typer.Option("--verify-hash", help="Verify the index's stored content hash against the MAF file (reads the whole file).")] = False, + progress_interval: OPTIONS.ProgressIntervalOption = COMMON.PROGRESS_INTERVAL_SECONDS, ) -> None: args = SimpleNamespace( maf_file=maf_file, @@ -2334,6 +2566,7 @@ def fetch_command( verbose=verbose, profile=profile, verify_hash=verify_hash, + progress_interval=progress_interval, ) cmdline = "mafutils fetch" if len(sys.argv) > 2: diff --git a/mafutils/gc.py b/mafutils/gc.py index 0a648c4..a0e97dd 100644 --- a/mafutils/gc.py +++ b/mafutils/gc.py @@ -51,6 +51,7 @@ from mafutils.lib import common as COMMON from mafutils.lib import loginit as LOGINIT +from mafutils.lib import options as OPTIONS ############################################################################# @@ -154,13 +155,26 @@ def chunker(iterable, size): ############################################################################# -def runSequentialGC(maf_file, maf_compression, LOG): +def runSequentialGC(maf_file, maf_compression, LOG, progress_interval=COMMON.PROGRESS_INTERVAL_SECONDS): counts = defaultdict(newBaseCounts) + # No block total without the index, so progress is position in the file + # on disk. Checked on "a" lines only, every PROGRESS_CHECK_BLOCKS blocks, + # so the s lines -- nearly all of the loop -- pay nothing for it; with + # --progress-interval 0 the "a" lines skip it too. + LOG.info(COMMON.progressStartupLine(progress_interval)) + progress = COMMON.ByteProgress(LOG, os.path.getsize(maf_file), interval=progress_interval) + report_progress = progress.enabled + n_blocks = 0 + with COMMON.openMaf(maf_file, maf_compression, "rt") as maf_fp: for line in maf_fp: if line.startswith("s"): updateCountsFromSeqLine(line, counts) + elif report_progress and line.startswith("a"): + n_blocks += 1 + if n_blocks % COMMON.PROGRESS_CHECK_BLOCKS == 0: + progress.updateTo(COMMON.diskPosition(maf_fp)) return counts @@ -196,9 +210,13 @@ def mergeCountsInto(target, counts): m["T"] += rec["T"] -def runParallelGC(maf_file, maf_compression, index_file, processes, chunk_size, LOG, verify_hash=False): +def runParallelGC(maf_file, maf_compression, index_file, processes, chunk_size, LOG, verify_hash=False, + progress_interval=COMMON.PROGRESS_INTERVAL_SECONDS): LOG.info(f"Streaming index file: {index_file}") - COMMON.validateIndexHeader(COMMON.readIndexHeader(index_file), maf_file, maf_compression, LOG, strict=verify_hash) + index_header = COMMON.readIndexHeader(index_file) + COMMON.validateIndexHeader(index_header, maf_file, maf_compression, LOG, strict=verify_hash) + header_blocks = COMMON.headerBlockCount(index_header, LOG, progress=progress_interval > 0) + LOG.info(COMMON.progressStartupLine(progress_interval)) LOG.info(f"Streaming index in chunks of {chunk_size} blocks across {processes} process(es).") @@ -213,6 +231,7 @@ def runParallelGC(maf_file, maf_compression, index_file, processes, chunk_size, total_blocks = 0 tasks_done = 0 chunk_iter = enumerate(chunker(iterIndexEntries(index_file, LOG), chunk_size), start=1) + progress = COMMON.BlockProgress(LOG, header_blocks, interval=progress_interval) with ProcessPoolExecutor(max_workers=processes) as executor: in_flight = {} @@ -224,20 +243,31 @@ def runParallelGC(maf_file, maf_compression, index_file, processes, chunk_size, for future in done: # Merge as results arrive, then drop them, so peak memory is # one accumulator plus the in-flight chunks. - total_blocks += in_flight.pop(future) + n_blocks = in_flight.pop(future) + total_blocks += n_blocks mergeCountsInto(merged, future.result()) tasks_done += 1 + progress.update(n_blocks) for future in in_flight: total_blocks += in_flight[future] mergeCountsInto(merged, future.result()) tasks_done += 1 + progress.update(in_flight[future]) if not total_blocks: LOG.error("No valid index entries found.") sys.exit(1) LOG.info(f"Processed {total_blocks} indexed blocks in {tasks_done} task(s).") + if header_blocks is not None and total_blocks != header_blocks: + # Same check as stats: both come from one `mafutils index` run, so they + # only disagree if the index was truncated or edited. + LOG.warning( + f"Processed {total_blocks:,} blocks but the index header records " + f"{header_blocks:,}. The index rows and header disagree, so this GC content may " + f"not cover the whole MAF; rebuild the index with `mafutils index`." + ) return merged @@ -304,8 +334,9 @@ def run_gc(args, cmdline="mafutils gc"): LOG.info(f"mafutils gc called as: {cmdline}") LOG.info("-" * 40) + arg_width = max(len(arg) for arg in vars(args)) for arg, value in vars(args).items(): - LOG.info(f"{arg:16} : {value}") + LOG.info(f"{arg:{arg_width}} : {value}") LOG.info("-" * 40) if args.processes < 1: @@ -347,11 +378,11 @@ def run_gc(args, cmdline="mafutils gc"): LOG.info(f"Running in PARALLEL mode ({maf_compression} MAF + index + processes > 1)") counts = runParallelGC( args.maf_file, maf_compression, index_file, args.processes, args.chunk_size, LOG, - verify_hash=args.verify_hash, + verify_hash=args.verify_hash, progress_interval=args.progress_interval, ) else: LOG.info("Running in SEQUENTIAL mode") - counts = runSequentialGC(args.maf_file, maf_compression, LOG) + counts = runSequentialGC(args.maf_file, maf_compression, LOG, progress_interval=args.progress_interval) species_gc = computeSpeciesGC(counts, LOG) if not species_gc: @@ -384,6 +415,7 @@ def gc_command( chunk_size: Annotated[int, typer.Option("--chunk-size", help="Blocks per worker task in parallel mode (default: 5000)")] = 5000, verify_hash: Annotated[bool, typer.Option("--verify-hash", help="Verify the index's stored content hash against the MAF file (reads the whole file; only checked when an index is actually used).")] = False, log_level: Annotated[LogLevel, typer.Option("--log-level", help="Logging level (default: INFO)")] = LogLevel.INFO, + progress_interval: OPTIONS.ProgressIntervalOption = COMMON.PROGRESS_INTERVAL_SECONDS, ) -> None: args = SimpleNamespace( maf_file=maf_file, @@ -393,6 +425,7 @@ def gc_command( chunk_size=chunk_size, verify_hash=verify_hash, log_level=log_level.value, + progress_interval=progress_interval, ) cmdline = "mafutils gc" if len(sys.argv) > 2: diff --git a/mafutils/index.py b/mafutils/index.py index 7963b00..d92962d 100644 --- a/mafutils/index.py +++ b/mafutils/index.py @@ -6,6 +6,7 @@ ############################################################################# import hashlib +import logging import os import shutil import tempfile @@ -14,6 +15,8 @@ import typer from mafutils.lib import common as COMMON +from mafutils.lib import loginit as LOGINIT +from mafutils.lib import options as OPTIONS def process_maf_block(block): @@ -35,12 +38,20 @@ def process_maf_block(block): return row, int(ref_seq[3]), seq_len, num_seqs -def run_index(maf_file, block_index_path, scaffold_index_path): +def run_index(maf_file, block_index_path, scaffold_index_path, LOG, progress_interval=COMMON.PROGRESS_INTERVAL_SECONDS): maf_compression = COMMON.detectCompression(maf_file) size = os.path.getsize(maf_file) mtime = os.path.getmtime(maf_file) hash_obj = hashlib.md5() + LOG.info( + f"Indexing {maf_file} ({COMMON.humanBytes(size)}, {maf_compression}). " + f"{COMMON.progressStartupLine(progress_interval)}" + ) + # Bytes, not blocks: the block total is what this run is building. + progress = COMMON.ByteProgress(LOG, size, interval=progress_interval) + report_progress = progress.enabled # a local: tested once per block + out_dir = os.path.dirname(block_index_path) or "." with tempfile.TemporaryDirectory(prefix="maf_index_tmp_", dir=out_dir) as tmp_dir: block_tmp_path = os.path.join(tmp_dir, "block.idx") @@ -76,6 +87,11 @@ def run_index(maf_file, block_index_path, scaffold_index_path): ref_scaffold = block_info[0] n_blocks += 1 + # Hot path: only every PROGRESS_CHECK_BLOCKS blocks does this + # cost more than one modulo (a clock read and an lseek), and + # with --progress-interval 0 not even the modulo runs. + if report_progress and n_blocks % COMMON.PROGRESS_CHECK_BLOCKS == 0: + progress.updateTo(COMMON.diskPosition(maf_stream)) total_ref_bases += ref_len total_aln_cols += aln_cols total_seq_lines += num_seqs @@ -124,14 +140,27 @@ def run_index(maf_file, block_index_path, scaffold_index_path): with open(scaffold_tmp_path, "r") as tmp_fp: shutil.copyfileobj(tmp_fp, scaffold_stream) + LOG.info( + f"Indexed {n_blocks:,} blocks on {len(distinct_scaffolds):,} reference scaffold(s): " + f"{block_index_path}, {scaffold_index_path}" + ) + def index_command( maf_file: Annotated[str, typer.Argument(help="Input MAF file (.maf, .maf.gz, or bgzip-compressed .maf)")], block_index: Annotated[Optional[str], typer.Argument(help="Output block index path (default: .block.idx)")] = None, scaffold_index: Annotated[Optional[str], typer.Argument(help="Output scaffold index path (default: .scaffold.idx)")] = None, + progress_interval: OPTIONS.ProgressIntervalOption = COMMON.PROGRESS_INTERVAL_SECONDS, ) -> None: if block_index is None: block_index = COMMON.deriveBlockIndexPath(maf_file) if scaffold_index is None: scaffold_index = COMMON.deriveScaffoldIndexPath(maf_file) - run_index(maf_file, block_index, scaffold_index) + # Screen only, like validate: a log file next to the indexes would be a + # new, unasked-for output. + logger_name = LOGINIT.configureLogging( + log_level="INFO", + log_verbosity="SCREEN", + logger_name="maf_index_logger", + ) + run_index(maf_file, block_index, scaffold_index, logging.getLogger(logger_name), progress_interval=progress_interval) diff --git a/mafutils/info.py b/mafutils/info.py index 8b2bc2d..3df5694 100644 --- a/mafutils/info.py +++ b/mafutils/info.py @@ -43,17 +43,6 @@ ############################################################################# -def humanBytes(n): - """ - Decimal units, so "GB" means 10^9 bytes and the printed number matches the - `size=` field in the index header rather than being 7% smaller. - """ - for unit in ("B", "KB", "MB", "GB", "TB"): - if n < 1000 or unit == "TB": - return f"{n:.2f} {unit}" if unit != "B" else f"{n} B" - n /= 1000 - - def readAggregates(header): """ Pulls the aggregate counts out of a parsed index header, or None if this @@ -161,7 +150,7 @@ def run_info(args, LOG): LOG.warning( "Index predates the aggregate header fields; computing block stats by " "scanning the index (%s). Rebuild with `mafutils index` to make this instant.", - humanBytes(index_size), + COMMON.humanBytes(index_size), ) aggregates = scanIndexAggregates(index_file) @@ -187,7 +176,7 @@ def run_info(args, LOG): size = os.path.getsize(args.maf_file) out = [] out.append(f"{'MAF':<15}: {os.path.basename(args.maf_file)}") - out.append(f"{'size':<15}: {humanBytes(size)} ({maf_compression})") + out.append(f"{'size':<15}: {COMMON.humanBytes(size)} ({maf_compression})") if not have_index: out.append(f"{'index':<15}: not found at {index_file}") diff --git a/mafutils/lib/common.py b/mafutils/lib/common.py index 8d86cf5..8dec4c2 100644 --- a/mafutils/lib/common.py +++ b/mafutils/lib/common.py @@ -9,6 +9,7 @@ import io import os import sys +import time from mafutils.lib import bgzf @@ -114,6 +115,11 @@ def seekable(self): def tell(self): return self._raw.tell() + def fileno(self): + # For diskPosition() progress reporting only; all reads still go + # through readinto() below, so nothing bypasses the hash. + return self._raw.fileno() + def seek(self, offset, whence=0): # bgzf.BgzfReader issues "seeks" as part of purely sequential reading # (e.g. re-asserting the current position when loading each new @@ -198,7 +204,8 @@ def computeFileHash(maf_file, algo="md5"): COPY_CHUNK_SIZE = 8 * 1024 * 1024 -def copyMafRangeToStream(handle, compression, offset_start, offset_end, out_stream, chunk_size=COPY_CHUNK_SIZE): +def copyMafRangeToStream(handle, compression, offset_start, offset_end, out_stream, chunk_size=COPY_CHUNK_SIZE, + on_progress=None): """ Copies the byte range [offset_start, offset_end) from an open binary-mode MAF handle straight into out_stream, in bounded chunks. Returns the number @@ -215,6 +222,11 @@ def copyMafRangeToStream(handle, compression, offset_start, offset_end, out_stre Memory here is O(chunk_size) regardless of range size, for every compression type -- including bgzip, whose line-by-line path must also write as it goes rather than accumulating lines to join at the end. + + on_progress, if given, is called with the cumulative bytes written after + each chunk -- once per chunk_size bytes, so free at any range size. The + bgzip path never calls it: that loop runs once per LINE, and is kept free + of any per-line work. """ handle.seek(offset_start) written = 0 @@ -239,6 +251,8 @@ def copyMafRangeToStream(handle, compression, offset_start, offset_end, out_stre out_stream.write(chunk) written += len(chunk) remaining -= len(chunk) + if on_progress is not None: + on_progress(written) return written @@ -563,3 +577,168 @@ def iterMafBlocks(stream, binary=False): yield block_lines, block_start, pos if binary else stream.tell() ############################################################################# + +# Seconds between progress lines: the default for --progress-interval, and also +# its minimum (0 turns progress off; see lib/options.py). A plain periodic log +# line rather than a redrawn progress bar: the log goes to a file (and to SLURM +# output) as well as the screen, where a redrawn bar turns into garbage. +PROGRESS_INTERVAL_SECONDS = 60 + +# Loops that see every block (index, sequential gc) only check progress every +# this many blocks: a clock read plus an lseek is cheap, but not next to a +# per-block hot path run 233M times on a whole-genome MAF. +PROGRESS_CHECK_BLOCKS = 10_000 + + +def humanBytes(n): + """ + Decimal units, so "GB" means 10^9 bytes and the printed number matches the + `size=` field in the index header rather than being 7% smaller. + """ + for unit in ("B", "KB", "MB", "GB", "TB"): + if n < 1000 or unit == "TB": + return f"{n:.2f} {unit}" if unit != "B" else f"{n} B" + n /= 1000 + + +def diskPosition(stream): + """ + How far into the file on disk this process has read, for any stream from + openMaf() or openMafHashing(). An lseek on the shared file descriptor, so + for gz/bgzip it is COMPRESSED bytes -- comparable with the file size, which + the stream's own positions (decompressed offsets for gz, virtual offsets + for bgzip) are not. Runs ahead of parsing by at most one read buffer. + """ + return os.lseek(stream.fileno(), 0, os.SEEK_CUR) + + +def headerBlockCount(header, LOG, progress=True): + """ + The block count recorded in an index header (indexes built by mafutils + >= 0.7), used as the progress total and for the end-of-run check against + the blocks actually processed -- or None. When progress is on, says why it + will lack a percentage, rather than leaving it to look broken. + """ + if header and "blocks" in header: + try: + return int(header["blocks"]) + except ValueError: + consequence = "; progress will be reported without a percentage or time estimate" if progress else "" + LOG.warning(f"Index header has a non-integer block count ({header['blocks']!r}){consequence}.") + return None + if progress: + LOG.info( + "Index header has no block count (built before mafutils 0.7), so progress will be " + "reported without a percentage or time estimate. Rebuild the index with " + "`mafutils index` to get them." + ) + return None + + +def progressStartupLine(interval): + """The one line each command logs at startup about --progress-interval.""" + if interval > 0: + return f"Progress is logged every {interval}s." + return "Progress reporting is off (--progress-interval 0)." + + +def formatDuration(seconds): + seconds = int(round(seconds)) + hours, rem = divmod(seconds, 3600) + minutes, secs = divmod(rem, 60) + if hours: + return f"{hours}h {minutes:02d}m" + if minutes: + return f"{minutes}m {secs:02d}s" + return f"{secs}s" + + +class BlockProgress: + """ + Throttled "N / total blocks" log lines for a long block-processing loop. + + update() is called from the main process as each chunk of blocks finishes, + so its cost is one clock read per chunk and workers are never touched. It + logs at most once per `interval` seconds. + + `total` is the block count from the index header (indexes built by + mafutils >= 0.7), or None for an older index -- lines then report blocks + and rate without a percentage or time estimate. Blocks are a rough proxy + for time (their size varies across a genome), so the estimate is rough. + + `unit` names what is counted ("regions" for fetch). A caller that waits on + worker futures passes dueIn() as the wait timeout, so it wakes exactly when + a line is due instead of polling. + + interval=0 (--progress-interval 0) turns it off: update/updateTo return at + once without reading the clock, and dueIn() is None, so a wait() given it + blocks exactly as it did before progress existed. Callers with per-block + hot loops also test `enabled` themselves, to skip their own bookkeeping. + """ + + def __init__(self, LOG, total, interval=None, clock=time.monotonic, unit="blocks"): + self.LOG = LOG + self.total = total + self.unit = unit + # Read at construction, not bound at definition, so a test can lower + # the module constant to see real lines without waiting a minute. + self.interval = PROGRESS_INTERVAL_SECONDS if interval is None else interval + self.enabled = self.interval > 0 + self.clock = clock + self.done = 0 + self.start = clock() + self.last_logged = self.start + + def update(self, n_blocks): + if self.enabled: + self.updateTo(self.done + n_blocks) + + def updateTo(self, done): + if not self.enabled: + return + self.done = done + now = self.clock() + if now - self.last_logged >= self.interval: + self.last_logged = now + self.LOG.info(self.message(now)) + + def dueIn(self): + if not self.enabled: + return None + return max(0.0, self.last_logged + self.interval - self.clock()) + + def describe(self, done, total): + return f"{done:,} / {total:,} {self.unit}" if total else f"{done:,} {self.unit}" + + def describeRate(self, rate): + return f"{rate:,.0f} {self.unit}/s" + + def message(self, now): + elapsed = now - self.start + rate = self.done / elapsed if elapsed > 0 else 0.0 + msg = f"Progress: {self.describe(self.done, self.total)}" + if self.total: + msg += f" ({100 * self.done / self.total:.1f}%)" + msg += f", {self.describeRate(rate)}, {formatDuration(elapsed)} elapsed" + if self.total and rate > 0 and self.done < self.total: + msg += f", ~{formatDuration((self.total - self.done) / rate)} left" + return msg + + +class ByteProgress(BlockProgress): + """ + BlockProgress measured in bytes of the file on disk, for a single + sequential pass over the MAF that has no block total to measure against + (`mafutils index` is what builds it; sequential `gc` doesn't read the + index). Fed with updateTo(diskPosition(stream)), so for gz/bgzip it is + compressed bytes against the compressed file size. Bytes track time more + evenly than blocks do, since parsing cost scales with bytes. + """ + + def describe(self, done, total): + return f"{humanBytes(done)} / {humanBytes(total)}" if total else humanBytes(done) + + def describeRate(self, rate): + return f"{humanBytes(int(rate))}/s" + +############################################################################# diff --git a/mafutils/lib/options.py b/mafutils/lib/options.py new file mode 100644 index 0000000..08800e1 --- /dev/null +++ b/mafutils/lib/options.py @@ -0,0 +1,42 @@ +############################################################################# +# Command-line options shared by several mafutils subcommands, defined once so +# their name, default, help and validation cannot drift apart between them. +############################################################################# + +from typing import Annotated + +import typer + +from mafutils.lib import common as COMMON + + +def validateProgressInterval(value: int) -> int: + """ + 0 (off) or at least the default. Anything else is a usage error rather than + being clamped -- a silently changed value would hide the mistake. The + minimum applies at the command line only; the run_* functions accept any + interval, so tests can use short ones. + """ + if value == 0 or value >= COMMON.PROGRESS_INTERVAL_SECONDS: + return value + raise typer.BadParameter( + f"must be 0 (off) or at least {COMMON.PROGRESS_INTERVAL_SECONDS} seconds, not {value}." + ) + + +# Per command (index, stats, gc, fetch) rather than one global option: a global +# option must come before the subcommand (`mafutils --progress-interval 300 +# stats ...`), so the natural `mafutils stats ... --progress-interval 300` would +# fail, and it would not appear in each command's --help. +ProgressIntervalOption = Annotated[ + int, + typer.Option( + "--progress-interval", + metavar="SECONDS", + callback=validateProgressInterval, + help=( + f"Seconds between progress lines (default and minimum: {COMMON.PROGRESS_INTERVAL_SECONDS}). " + f"0 turns progress reporting off." + ), + ), +] diff --git a/mafutils/stats.py b/mafutils/stats.py index fd68906..77bdec6 100644 --- a/mafutils/stats.py +++ b/mafutils/stats.py @@ -13,6 +13,8 @@ Outputs: .overall.tsv : one metric per line .species.tsv : one row per species + .species_scaffolds.tsv : one row per (species, source scaffold) + present in the MAF -- not the whole assembly .block.tsv : one row per block (optional; default enabled) .dashboard.html : self-contained summary dashboard (optional) @@ -42,6 +44,7 @@ from mafutils.lib import common as COMMON from mafutils.lib import loginit as LOGINIT +from mafutils.lib import options as OPTIONS ############################################################################# @@ -140,8 +143,44 @@ def speciesFromSrc(src): return src +class SrcSizeConflictError(ValueError): + """Raised when one species.scaffold reports two different srcSize values.""" + + +def reconcileSrcSize(key, known, incoming): + """ + Combines two srcSize observations for the same (species, scaffold). + + srcSize is the scaffold's full length in its assembly, so every s line for + one species.scaffold must agree. Disagreement means an internally + inconsistent MAF and is raised rather than resolved: picking either value + would be a silent fallback, and "whichever was seen first" would also make + the output depend on how blocks were split across workers. (Measured on + real data: zero conflicts across 2,647 scaffolds / 500,000 blocks.) + + None means "this line's srcSize did not parse" and never conflicts. + """ + if known is None: + return incoming + if incoming is None or incoming == known: + return known + species, scaffold = key + raise SrcSizeConflictError( + f"{species}.{scaffold} has conflicting srcSize values ({known} and {incoming}); " + f"the MAF is internally inconsistent, so its scaffold lengths cannot be reported." + ) + + def parseBlock(block_text): species_stats = defaultdict(lambda: {"copies": 0, "gaps": 0, "nongaps": 0}) + # (species, scaffold) -> {"src_size", "aligned_bases"} for this block. + # Collected here because this loop already splits every s line, so the + # scaffold name and srcSize (fields[1], fields[5]) are already in hand -- + # unlike `mafutils index`, which deliberately never parses non-reference + # lines and would pay ~2x to collect them. Keyed by species AND scaffold: + # same-named scaffolds in different species are different sequences. + scaffold_stats = {} + scaffold_size_parse_errors = 0 all_gap_species = set() present_species = set() n_seq_lines = 0 @@ -158,7 +197,9 @@ def parseBlock(block_text): if len(fields) < 7: parse_errors += 1 continue - sp = speciesFromSrc(fields[1]) + # One partition yields both halves; the species half is identical to + # speciesFromSrc(), so existing per-species statistics are unchanged. + sp, _dot, scaffold = fields[1].partition(".") seq = fields[6] gaps = seq.count("-") nongaps = len(seq) - gaps @@ -171,6 +212,23 @@ def parseBlock(block_text): rec["gaps"] += gaps rec["nongaps"] += nongaps + # A non-integer srcSize affects ONLY the scaffold table: this line has + # already been counted toward every existing statistic above, exactly + # as before (which never read fields[5]). Its scaffold is still + # recorded, without a size, and the failure is counted -- not hidden. + try: + src_size = int(fields[5]) + except ValueError: + src_size = None + scaffold_size_parse_errors += 1 + key = (sp, scaffold) + srec = scaffold_stats.get(key) + if srec is None: + scaffold_stats[key] = {"src_size": src_size, "aligned_bases": nongaps} + else: + srec["aligned_bases"] += nongaps + srec["src_size"] = reconcileSrcSize(key, srec["src_size"], src_size) + for sp, rec in species_stats.items(): if rec["nongaps"] == 0: all_gap_species.add(sp) @@ -211,6 +269,8 @@ def parseBlock(block_text): return { "species_stats": species_stats, + "scaffold_stats": scaffold_stats, + "scaffold_size_parse_errors": scaffold_size_parse_errors, "present_species": present_species, "all_gap_species": all_gap_species, "n_seq_lines": n_seq_lines, @@ -258,6 +318,7 @@ def workerTask( "total_gap_chars": 0, "total_cells": 0, "total_parse_errors": 0, + "scaffold_size_parse_errors": 0, "blocks_with_zero_parsed_seq_lines_but_index_nonzero": 0, "all_gap_sites": 0, "invariant_sites": 0, @@ -278,6 +339,9 @@ def workerTask( } ) observed_species = set() + # (species, scaffold) -> {"src_size", "blocks", "aligned_bases"} for this + # chunk; merged into the run-wide accumulator as each chunk completes. + scaffolds = {} block_tmp = None block_tmp_path = None @@ -325,6 +389,20 @@ def workerTask( observed_species.update(parsed["present_species"]) + overall["scaffold_size_parse_errors"] += parsed["scaffold_size_parse_errors"] + for key, srec in parsed["scaffold_stats"].items(): + crec = scaffolds.get(key) + if crec is None: + scaffolds[key] = { + "src_size": srec["src_size"], + "blocks": 1, + "aligned_bases": srec["aligned_bases"], + } + else: + crec["blocks"] += 1 + crec["aligned_bases"] += srec["aligned_bases"] + crec["src_size"] = reconcileSrcSize(key, crec["src_size"], srec["src_size"]) + for sp, rec in parsed["species_stats"].items(): s = species[sp] s["blocks_present"] += 1 @@ -373,6 +451,7 @@ def workerTask( "overall": overall, "species": dict(species), "observed_species": sorted(observed_species), + "scaffolds": scaffolds, "block_tmp_path": block_tmp_path, } @@ -401,6 +480,10 @@ def newStatsAccumulator(): ), "observed_species": set(), "block_tmp_paths": [], + # (species, scaffold) -> {"src_size", "blocks", "aligned_bases"}. + # Bounded by the number of distinct scaffolds, not by blocks: 1,587 + # pairs measured at ~0.1MB over 200,000 real blocks. + "scaffolds": {}, } @@ -419,6 +502,16 @@ def mergeStatsInto(acc, res): for k, v in rec.items(): s[k] += v acc["observed_species"].update(res["observed_species"]) + scaffolds = acc["scaffolds"] + for key, crec in res["scaffolds"].items(): + arec = scaffolds.get(key) + if arec is None: + scaffolds[key] = dict(crec) + else: + arec["blocks"] += crec["blocks"] + arec["aligned_bases"] += crec["aligned_bases"] + # Raises on a genuine conflict; see reconcileSrcSize. + arec["src_size"] = reconcileSrcSize(key, arec["src_size"], crec["src_size"]) if res["block_tmp_path"]: # Keyed by task_id so the per-chunk block rows can be concatenated in # index order. Results arrive out of order in the parallel path, and @@ -453,6 +546,10 @@ def writeOverall(out_path, overall, n_species_total): "blocks_with_zero_parsed_seq_lines_but_index_nonzero", int(overall["blocks_with_zero_parsed_seq_lines_but_index_nonzero"]), ), + # s lines whose srcSize did not parse as an integer. Affects only the + # species_scaffolds table (those scaffolds are reported without a + # size); every other statistic counted the line as before. + ("scaffold_size_parse_errors", int(overall["scaffold_size_parse_errors"])), ] with open(out_path, "w", encoding="utf-8") as fp: @@ -464,7 +561,28 @@ def writeOverall(out_path, overall, n_species_total): fp.write(f"{key}\t{val}\n") -def writeSpecies(out_path, species, total_blocks): +def speciesScaffoldRollup(scaffolds): + """ + species -> (scaffolds_in_maf, scaffold_bases_in_maf). + + Both are limited to scaffolds PRESENT IN THE MAF, i.e. with at least one + aligned s line. A scaffold that never aligns to the reference never + appears, so these are lower bounds on the assembly's scaffold count and + size -- while also over-counting relative to aligned bases, since a + barely-aligned scaffold contributes its full srcSize. Scaffolds whose + srcSize failed to parse are counted but contribute no bases. + """ + rollup = defaultdict(lambda: [0, 0]) + for (species, _scaffold), rec in scaffolds.items(): + r = rollup[species] + r[0] += 1 + if rec["src_size"] is not None: + r[1] += rec["src_size"] + return rollup + + +def writeSpecies(out_path, species, total_blocks, scaffolds=None): + rollup = speciesScaffoldRollup(scaffolds or {}) header = [ "species", "blocks_present", @@ -479,6 +597,11 @@ def writeSpecies(out_path, species, total_blocks): "avg_gaps_per_present_block", "avg_gap_fraction_per_present_block", "all_gap_blocks", + # Appended last so positional parsers of species.tsv keep working. + # "_in_maf" because these are NOT assembly properties -- see + # speciesScaffoldRollup. + "scaffolds_in_maf", + "scaffold_bases_in_maf", ] with open(out_path, "w", encoding="utf-8") as fp: fp.write("\t".join(header) + "\n") @@ -503,12 +626,53 @@ def writeSpecies(out_path, species, total_blocks): f"{(rec['gaps_total'] / blocks_present) if blocks_present else 0.0:.8f}", f"{(rec['gap_pct_block_sum'] / blocks_present) if blocks_present else 0.0:.8f}", str(int(rec["all_gap_blocks"])), + str(rollup[sp][0]) if sp in rollup else "0", + str(rollup[sp][1]) if sp in rollup else "0", ] ) + "\n" ) +SPECIES_SCAFFOLD_HEADERS = [ + "species", "scaffold", "src_size", "blocks", "aligned_bases", "aligned_fraction", +] + + +def writeSpeciesScaffolds(out_path, scaffolds): + """ + One row per (species, scaffold) present in the MAF, sorted by species then + src_size descending (largest scaffolds first). + + src_size is the scaffold's true full length in its assembly (MAF srcSize), + but only scaffolds with at least one aligned s line appear at all -- the + table is complete over the MAF, not over the assembly. + + aligned_fraction is deliberately not clamped: a region aligned into several + blocks (duplicated/paralogous alignments) can push it above 1.0, which is + information rather than an error. It is "." when src_size did not parse or + is zero. + """ + def sort_key(item): + (species, scaffold), rec = item + size = rec["src_size"] if rec["src_size"] is not None else -1 + return (species, -size, scaffold) + + with open(out_path, "w", encoding="utf-8") as fp: + fp.write("\t".join(SPECIES_SCAFFOLD_HEADERS) + "\n") + for (species, scaffold), rec in sorted(scaffolds.items(), key=sort_key): + size = rec["src_size"] + fraction = f"{rec['aligned_bases'] / size:.8f}" if size else "." + fp.write("\t".join([ + species, + scaffold or ".", + "." if size is None else str(size), + str(rec["blocks"]), + str(rec["aligned_bases"]), + fraction, + ]) + "\n") + + def writeBlockTable(out_path, block_tmp_paths, n_species_total, expected_species): expected_species = set(expected_species) if expected_species else None header = [ @@ -646,6 +810,8 @@ def writeDashboard( cmdline, generated_at, LOG, + scaffolds=None, + top_scaffolds=10, ): try: import matplotlib @@ -668,10 +834,14 @@ def writeDashboard( parsimony_informative_sites = int(overall["parsimony_informative_sites"]) all_gap_sites = int(overall["all_gap_sites"]) + scaffolds = scaffolds or {} + rollup = speciesScaffoldRollup(scaffolds) + rows = [] for sp, rec in species.items(): blocks_present = int(rec["blocks_present"]) total_sp_cells = rec["gaps_total"] + rec["nongaps_total"] + n_scaf, scaf_bases = rollup[sp] if sp in rollup else (0, 0) rows.append( { "species": sp, @@ -679,6 +849,8 @@ def writeDashboard( "gap_fraction_cells": (rec["gaps_total"] / total_sp_cells) if total_sp_cells else 0.0, "duplicated_blocks": int(rec["duplicated_blocks"]), "blocks_present": blocks_present, + "scaffolds_in_maf": n_scaf, + "scaffold_bases_in_maf": scaf_bases, } ) @@ -857,6 +1029,44 @@ def makeBarFigure(data_rows, x_key, title, color): ) species_summary_html = "".join(f"
  • {html.escape(item)}
    • " for item in species_summary_items) + + # Largest scaffolds per species: one collapsed
    per species, so a + # 241-way MAF stays scannable. Each summary line states "showing N of M" + # up front -- a trailing "... and N more" is easy to miss (as found with + # the species list in `mafutils info`). The TSV is never capped. + top_scaffolds = max(1, int(top_scaffolds)) + by_species = defaultdict(list) + for (sp, scaffold), rec in scaffolds.items(): + by_species[sp].append((scaffold, rec)) + scaffold_detail_parts = [] + for sp in sorted(by_species): + items = sorted( + by_species[sp], + key=lambda kv: (-(kv[1]["src_size"] if kv[1]["src_size"] is not None else -1), kv[0]), + ) + shown = items[:top_scaffolds] + body = "".join( + "" + f"{html.escape(scaffold or '.')}" + f"{'.' if rec['src_size'] is None else format(rec['src_size'], ',')}" + f"{rec['blocks']:,}" + f"{rec['aligned_bases']:,}" + f"{(format(rec['aligned_bases'] / rec['src_size'], '.4f')) if rec['src_size'] else '.'}" + "" + for scaffold, rec in shown + ) + scaffold_detail_parts.append( + "
    " + f"{html.escape(sp)} — showing {len(shown):,} of {len(items):,} scaffolds" + f"{' (largest by length)' if len(items) > len(shown) else ''}" + "
    " + "" + "" + "" + f"{body}" + "
    ScaffoldLength (src_size)BlocksAligned basesAligned fraction
    " + ) + scaffold_details_html = "\n".join(scaffold_detail_parts) or "

    No scaffolds recorded.

    " species_table_html = "\n".join( [ ( @@ -867,6 +1077,8 @@ def makeBarFigure(data_rows, x_key, title, color): f"{r['presence_fraction']:.4f}" f"{r['gap_fraction_cells']:.4f}" f"{r['duplicated_blocks']:,}" + f"{r['scaffolds_in_maf']:,}" + f"{r['scaffold_bases_in_maf']:,}" "" ) for r in species_table_rows @@ -1089,6 +1301,21 @@ def makeBarFigure(data_rows, x_key, title, color): position: sticky; top: 0; }} +details {{ + margin: 6px 0; + border: 1px solid #dfe6ee; + border-radius: 6px; + padding: 4px 10px; +}} +summary {{ + cursor: pointer; + padding: 4px 0; +}} +.note {{ + font-size: 12px; + color: #556270; + margin: 8px 0 4px; +}} .table-wrap {{ max-height: 320px; overflow: auto; @@ -1154,6 +1381,8 @@ def makeBarFigure(data_rows, x_key, title, color): Presence fraction Gap fraction (cells) Duplicated blocks + Scaffolds in MAF + Scaffold bases in MAF @@ -1161,6 +1390,18 @@ def makeBarFigure(data_rows, x_key, title, color): +

    Scaffold counts and bases cover only scaffolds present in the + MAF (at least one aligned sequence line), not the species' full assembly: + unaligned scaffolds never appear, so these are lower bounds on assembly size. + Each scaffold's length is its true full length (MAF srcSize).

    + +
    + Largest scaffolds per species +

    Scaffolds present in the MAF, largest first. Collapsed by + default; the full list for every species is in + <prefix>.species_scaffolds.tsv. Aligned fraction can exceed + 1.0 where a region aligns into several blocks.

    + {scaffold_details_html}
    {plot_section_html} {block_section_html} @@ -1216,8 +1457,9 @@ def run_stats(args, cmdline="mafutils stats"): LOG.info(f"mafutils stats called as: {cmdline}") LOG.info("-" * 40) + arg_width = max(len(arg) for arg in vars(args)) for arg, value in vars(args).items(): - LOG.info(f"{arg:22} : {value}") + LOG.info(f"{arg:{arg_width}} : {value}") LOG.info("-" * 40) if args.processes < 1: @@ -1248,7 +1490,11 @@ def run_stats(args, cmdline="mafutils stats"): LOG.info(f"Loaded {len(expected_species)} expected species for exact missing-species reporting.") LOG.info(f"Reading index file: {index_file}") - COMMON.validateIndexHeader(COMMON.readIndexHeader(index_file), args.maf_file, maf_compression, LOG, strict=args.verify_hash) + index_header = COMMON.readIndexHeader(index_file) + COMMON.validateIndexHeader(index_header, args.maf_file, maf_compression, LOG, strict=args.verify_hash) + + header_blocks = COMMON.headerBlockCount(index_header, LOG, progress=args.progress_interval > 0) + LOG.info(COMMON.progressStartupLine(args.progress_interval)) effective_processes = args.processes effective_chunk_size = args.chunk_size @@ -1276,66 +1522,90 @@ def run_stats(args, cmdline="mafutils stats"): tasks_done = 0 # The index is streamed, never materialized -- see iterIndexEntries. chunk_iter = enumerate(chunker(iterIndexEntries(index_file, LOG), effective_chunk_size), start=1) - - if effective_processes > 1: - # Bound in-flight work and merge results as they arrive, so peak - # memory is one accumulator plus the queued chunks rather than - # every chunk and every chunk's result. - max_in_flight = effective_processes * 2 - with ProcessPoolExecutor(max_workers=effective_processes) as executor: - in_flight = {} - for task_id, chunk in chunk_iter: - future = executor.submit( - workerTask, - task_id, - args.maf_file, - maf_compression, - chunk, - write_block_rows, - tmp_dir, - ) - in_flight[future] = len(chunk) - if len(in_flight) < max_in_flight: - continue - done, _pending = wait(list(in_flight.keys()), return_when=FIRST_COMPLETED) - for f in done: - total_blocks_seen += in_flight.pop(f) + progress = COMMON.BlockProgress(LOG, header_blocks, interval=args.progress_interval) + + try: + if effective_processes > 1: + # Bound in-flight work and merge results as they arrive, so peak + # memory is one accumulator plus the queued chunks rather than + # every chunk and every chunk's result. + max_in_flight = effective_processes * 2 + with ProcessPoolExecutor(max_workers=effective_processes) as executor: + in_flight = {} + for task_id, chunk in chunk_iter: + future = executor.submit( + workerTask, + task_id, + args.maf_file, + maf_compression, + chunk, + write_block_rows, + tmp_dir, + ) + in_flight[future] = len(chunk) + if len(in_flight) < max_in_flight: + continue + done, _pending = wait(list(in_flight.keys()), return_when=FIRST_COMPLETED) + for f in done: + n_blocks = in_flight.pop(f) + total_blocks_seen += n_blocks + mergeStatsInto(acc, f.result()) + tasks_done += 1 + progress.update(n_blocks) + for f in list(in_flight): + n_blocks = in_flight.pop(f) + total_blocks_seen += n_blocks mergeStatsInto(acc, f.result()) tasks_done += 1 - for f in list(in_flight): - total_blocks_seen += in_flight.pop(f) - mergeStatsInto(acc, f.result()) - else: - # Skip ProcessPoolExecutor entirely at processes=1 -- pool - # creation carries a large, mostly-fixed memory cost regardless - # of worker count (confirmed via mafutils gc's equivalent - # parallel/sequential split -- see DEVELOPMENT.md), not worth - # paying for a single worker that does the exact same work - # in-process anyway. workerTask is a plain function with no - # pool-initializer dependency, so calling it directly here is - # equivalent to submitting it to a 1-worker pool. - # - # One handle is opened here and reused for every chunk: entries - # are in ascending file order, so all seeks stay forward. This is - # what preserves the cheap gz path now that the index arrives in - # many chunks instead of one -- reopening per chunk would restart - # gzip decompression from byte 0 each time. - with COMMON.openMaf(args.maf_file, maf_compression, "rb") as shared_fp: - for task_id, chunk in chunk_iter: - total_blocks_seen += len(chunk) - tasks_done += 1 - mergeStatsInto( - acc, - workerTask( - task_id, args.maf_file, maf_compression, chunk, - write_block_rows, tmp_dir, maf_fp=shared_fp, - ), - ) + progress.update(n_blocks) + else: + # Skip ProcessPoolExecutor entirely at processes=1 -- pool + # creation carries a large, mostly-fixed memory cost regardless + # of worker count (confirmed via mafutils gc's equivalent + # parallel/sequential split -- see DEVELOPMENT.md), not worth + # paying for a single worker that does the exact same work + # in-process anyway. workerTask is a plain function with no + # pool-initializer dependency, so calling it directly here is + # equivalent to submitting it to a 1-worker pool. + # + # One handle is opened here and reused for every chunk: entries + # are in ascending file order, so all seeks stay forward. This is + # what preserves the cheap gz path now that the index arrives in + # many chunks instead of one -- reopening per chunk would restart + # gzip decompression from byte 0 each time. + with COMMON.openMaf(args.maf_file, maf_compression, "rb") as shared_fp: + for task_id, chunk in chunk_iter: + total_blocks_seen += len(chunk) + tasks_done += 1 + mergeStatsInto( + acc, + workerTask( + task_id, args.maf_file, maf_compression, chunk, + write_block_rows, tmp_dir, maf_fp=shared_fp, + ), + ) + progress.update(len(chunk)) + + except SrcSizeConflictError as exc: + # Raised by a worker (re-raised through future.result()) or by the + # merge. Reported like every other fatal condition here -- logged, + # so it also lands in the run's log file -- rather than as a traceback. + LOG.error(str(exc)) + sys.exit(1) if not total_blocks_seen: LOG.error("No valid index entries found.") sys.exit(1) LOG.info(f"Processed {total_blocks_seen} indexed blocks in {tasks_done} tasks.") + if header_blocks is not None and total_blocks_seen != header_blocks: + # Same run of `mafutils index` writes both, so they only disagree + # if the index was truncated or edited (or rows were skipped as + # malformed, each already warned about above). + LOG.warning( + f"Processed {total_blocks_seen:,} blocks but the index header records " + f"{header_blocks:,}. The index rows and header disagree, so these statistics may " + f"not cover the whole MAF; rebuild the index with `mafutils index`." + ) overall = acc["overall"] species = acc["species"] @@ -1376,13 +1646,31 @@ def run_stats(args, cmdline="mafutils stats"): overall_file = out_prefix + ".overall.tsv" species_file = out_prefix + ".species.tsv" + species_scaffolds_file = out_prefix + ".species_scaffolds.tsv" block_file = out_prefix + ".block.tsv" + scaffolds = acc["scaffolds"] LOG.info(f"Writing overall stats: {overall_file}") writeOverall(overall_file, overall, species_total) LOG.info(f"Writing species stats: {species_file}") - writeSpecies(species_file, species, int(overall["total_blocks"])) + writeSpecies(species_file, species, int(overall["total_blocks"]), scaffolds) + + LOG.info( + f"Writing per-species scaffolds: {species_scaffolds_file} " + f"({len(scaffolds):,} species/scaffold pairs present in the MAF)" + ) + writeSpeciesScaffolds(species_scaffolds_file, scaffolds) + + size_errors = int(overall["scaffold_size_parse_errors"]) + if size_errors: + LOG.warning( + "%d s line(s) had a srcSize that did not parse as an integer. A scaffold with no " + "parseable srcSize on any line is reported without a size in %s; all other " + "statistics are unaffected.", + size_errors, + species_scaffolds_file, + ) if args.html_dashboard: dashboard_file = out_prefix + ".dashboard.html" @@ -1398,6 +1686,8 @@ def run_stats(args, cmdline="mafutils stats"): cmdline, generated_at, LOG, + scaffolds=scaffolds, + top_scaffolds=args.dashboard_top_scaffolds, ) if not args.no_block_table: @@ -1423,8 +1713,10 @@ def stats_command( expected_species_file: Annotated[str, typer.Option("--expected-species-file", help="File with one species name per line for exact missing lists.")] = "", verify_hash: Annotated[bool, typer.Option("--verify-hash", help="Verify the index's stored content hash against the MAF file (reads the whole file).")] = False, log_level: Annotated[LogLevel, typer.Option("--log-level", help="Logging level (default: INFO)")] = LogLevel.INFO, + progress_interval: OPTIONS.ProgressIntervalOption = COMMON.PROGRESS_INTERVAL_SECONDS, html_dashboard: Annotated[bool, typer.Option("--html-dashboard", help="Write an HTML dashboard with summary metrics and species plots.")] = False, dashboard_top_species: Annotated[int, typer.Option("--dashboard-top-species", help="Number of species to show in each dashboard bar plot (default: 25).")] = 25, + dashboard_top_scaffolds: Annotated[int, typer.Option("--dashboard-top-scaffolds", help="Largest scaffolds listed per species in the dashboard (default: 10). The species_scaffolds.tsv table is never capped.")] = 10, dashboard_max_block_points: Annotated[int, typer.Option("--dashboard-max-block-points", help="Maximum block rows sampled for dashboard block plots (default: 50000).")] = 50000, ) -> None: args = SimpleNamespace( @@ -1438,8 +1730,10 @@ def stats_command( expected_species_file=expected_species_file, verify_hash=verify_hash, log_level=log_level.value, + progress_interval=progress_interval, html_dashboard=html_dashboard, dashboard_top_species=dashboard_top_species, + dashboard_top_scaffolds=dashboard_top_scaffolds, dashboard_max_block_points=dashboard_max_block_points, ) cmdline = "mafutils stats" diff --git a/tests/example-tabs.maf b/tests/example-tabs.maf index 3ea0c57..ae8c78f 100644 --- a/tests/example-tabs.maf +++ b/tests/example-tabs.maf @@ -37,4 +37,4 @@ s gorilla.chrX 99 4 - 800 AGTC a score=400 s human.chr4 30 5 + 1000 T--ACGT -s chimp.chr4 28 5 + 998 TGG-C-T +s chimp.chr4 28 5 + 1000 TGG-C-T diff --git a/tests/example-tabs.maf.block.idx b/tests/example-tabs.maf.block.idx index 0126974..b0c06c7 100644 --- a/tests/example-tabs.maf.block.idx +++ b/tests/example-tabs.maf.block.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example-tabs.maf compression=none size=777 mtime=1784145266.0 hash=md5:790612b936fb9a12ec1308268831db18 +# mafutils-index format=2 maf=example-tabs.maf compression=none size=778 mtime=1792526707.0 hash=md5:48230e6b5535b46c5e808c40e2dbe3f2 blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 chr1 0 10 10 35 2 16 101 chr1 18 8 8 33 2 101 182 chr1 32 8 8 33 2 182 263 @@ -8,4 +8,4 @@ chr3 6 4 4 28 2 411 482 chr4 0 5 5 29 2 482 557 chr4 6 4 4 28 2 557 628 chrX 50 4 4 28 2 628 700 -chr4 30 5 7 32 2 700 777 +chr4 30 5 7 32 2 700 778 diff --git a/tests/example-tabs.maf.scaffold.idx b/tests/example-tabs.maf.scaffold.idx index a8d6e61..d5e7762 100644 --- a/tests/example-tabs.maf.scaffold.idx +++ b/tests/example-tabs.maf.scaffold.idx @@ -1,7 +1,7 @@ -# mafutils-index format=2 maf=example-tabs.maf compression=none size=777 mtime=1784145266.0 hash=md5:790612b936fb9a12ec1308268831db18 +# mafutils-index format=2 maf=example-tabs.maf compression=none size=778 mtime=1792526707.0 hash=md5:48230e6b5535b46c5e808c40e2dbe3f2 blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 chr1 16 263 chr2 263 338 chr3 338 482 chr4 482 628 chrX 628 700 -chr4 700 777 +chr4 700 778 diff --git a/tests/example.maf b/tests/example.maf index 4c7e8c7..9481160 100644 --- a/tests/example.maf +++ b/tests/example.maf @@ -37,4 +37,4 @@ s gorilla.chrX 99 4 - 800 AGTC a score=400 s human.chr4 30 5 + 1000 T--ACGT -s chimp.chr4 28 5 + 998 TGG-C-T \ No newline at end of file +s chimp.chr4 28 5 + 1000 TGG-C-T \ No newline at end of file diff --git a/tests/example.maf.bgz b/tests/example.maf.bgz index 658898a5deb8aa45be4c1254efb09eadf589f692..c16374c0ca1a22662d44bdc1a41f2abd30a74584 100644 GIT binary patch literal 305 zcmb2|=3rp}f&Xj_PR>jWB8n&-dguR&giL_prGyAh&4gm%U4ax@9>y^zV{V#pN^u7myNj8UQWU0 z?}o=>R-Dydc(PUM&YWw-v+{(Qwx>$e8_rhx7|*;1v8p-a9mdMw6! zF!^0!g4v1o{jZb%zxr(2bN6eM{-621Z&tf~t6aS5(KO>7C5slwpOs2J@q5m~RcSGL s6<=@E${t>?_F`_*8SR&oSOr-3US(%r%FD*d%)lUz7AVpTKtU(~05U3nx&QzG literal 338 zcmb2|=3rp}f&Xj_PR>jWdJK#L4DbIdE9WMrDU_ua6=&w>*&1?9UFdOB0H{hHO*tE5 z?*!Yt!v;KV<@1}W4>tYsUb~R_V59hh30DoetQ;Z{EBs2d~^e zzCGOCuY~kylcP6*r;5=GTG!paGc9j&2zkn!_+fTH=!~;##C{d1@NB){JS)@w@)@O|TdN}LpM^fn3T0^8 z5|fh_V{5|T7!z)yH|tMa$*lmvYqusVrsLx?CmEtF~6uK!*#{Pnr;&Uw`CZAF`pn;85)jT4yIN=ajonp)@>5>Hy1-;b=2c-=L2?GECE02K; literal 285 zcmV+&0pk82iwFoF{#a@P|7Cb#ZE$R5E^T3E0F9D8Zo@DPg?FEVK*mCaep6}4AZXwP zG-Nf5v^EgOg^lL?p 0 and calls == [] + + +def _fetch(args): + return subprocess.run([sys.executable, "-m", "mafutils", "fetch"] + args, + check=False, cwd=REPO_ROOT, capture_output=True, text=True) + + +def test_fetch_block_mode_reports_progress_setup(tmp_path): + r = _fetch([MAF_FILE, os.path.join(TEST_DIR, "example.bed"), "-o", str(tmp_path), "-p", "2"]) + assert r.returncode == 0, r.stderr + assert "Progress is logged every 60s." in r.stdout + + +@pytest.mark.parametrize("fixture", ["example.maf", "example.maf.gz", "example.maf.bgz"]) +def test_fetch_scaffold_mode_reports_each_finished_scaffold(fixture, tmp_path): + bed = os.path.join(str(tmp_path), "s.bed") + with open(bed, "w") as fp: + fp.write("chr1\t0\t50\tchr1\nchr4\t0\t50\tchr4\n") + maf = os.path.join(TEST_DIR, fixture) + r = _fetch([maf, bed, "-i", maf + ".scaffold.idx", "-m", "scaffold", "-p", "2", + "-o", os.path.join(str(tmp_path), "out")]) + assert r.returncode == 0, r.stdout + r.stderr + assert "of 2 scaffolds)" in r.stdout + assert "Finished chr1" in r.stdout and "Finished chr4" in r.stdout + assert ("counts finished scaffolds" in r.stdout) == (fixture.endswith(".bgz")) + + +# --------------------------------------------------------------------------- +# --progress-interval: default and minimum 60, 0 = off, anything else rejected +# --------------------------------------------------------------------------- + + +class CountingClock(FakeClock): + def __init__(self): + super().__init__() + self.calls = 0 + + def __call__(self): + self.calls += 1 + return self.now + + +def test_interval_zero_turns_progress_off_without_reading_the_clock(log): + logger, messages = log + clock = CountingClock() + progress = COMMON.BlockProgress(logger, 1000, interval=0, clock=clock) + calls_at_construction = clock.calls + for _ in range(100): + clock.now += 3600 + progress.update(10) + progress.updateTo(500) + assert not progress.enabled + assert clock.calls == calls_at_construction # nothing ran + assert messages == [] and progress.done == 0 + assert progress.dueIn() is None # wait(timeout=None): plain blocking wait + + +def test_progress_startup_line(): + assert COMMON.progressStartupLine(60) == "Progress is logged every 60s." + assert COMMON.progressStartupLine(300) == "Progress is logged every 300s." + assert COMMON.progressStartupLine(0) == "Progress reporting is off (--progress-interval 0)." + + +def test_header_block_count_is_quiet_about_progress_when_it_is_off(log): + logger, messages = log + assert COMMON.headerBlockCount({"size": "777"}, logger, progress=False) is None + assert messages == [] + # a corrupt count is still worth a warning, just without the progress clause + assert COMMON.headerBlockCount({"blocks": "lots"}, logger, progress=False) is None + assert messages == ["Index header has a non-integer block count ('lots')."] + + +def _command_args(command, tmp_path): + out = str(tmp_path) + return { + "index": [MAF_FILE, os.path.join(out, "x.block.idx"), os.path.join(out, "x.scaffold.idx")], + "stats": [MAF_FILE, "-o", os.path.join(out, "s")], + "gc": [MAF_FILE, "-o", os.path.join(out, "g")], + "fetch": [MAF_FILE, os.path.join(TEST_DIR, "example.bed"), "-b", "id", "-o", os.path.join(out, "f")], + }[command] + + +def _mafutils(args): + return subprocess.run([sys.executable, "-m", "mafutils"] + args, + check=False, cwd=REPO_ROOT, capture_output=True, text=True) + + +COMMANDS = ["index", "stats", "gc", "fetch"] + + +@pytest.mark.parametrize("command", COMMANDS) +@pytest.mark.parametrize("extra,expected", [ + ([], "Progress is logged every 60s."), + (["--progress-interval", "120"], "Progress is logged every 120s."), + (["--progress-interval", "0"], "Progress reporting is off (--progress-interval 0)."), +]) +def test_progress_interval_values_accepted(command, extra, expected, tmp_path): + r = _mafutils([command] + _command_args(command, tmp_path) + extra) + assert r.returncode == 0, r.stdout + r.stderr + assert expected in r.stdout + assert "Progress:" not in r.stdout # runs this short never reach a line + + +@pytest.mark.parametrize("command", COMMANDS) +@pytest.mark.parametrize("value", ["30", "59", "-1"]) +def test_progress_interval_below_minimum_is_a_usage_error(command, value, tmp_path): + """Rejected, never clamped: a silently changed value would hide the mistake.""" + r = _mafutils([command] + _command_args(command, tmp_path) + ["--progress-interval", value]) + assert r.returncode == 2 + # typer draws the error in a box and wraps it; compare with the box and + # line breaks collapsed + flat = " ".join((r.stdout + r.stderr).replace("\u2502", " ").split()) + assert f"must be 0 (off) or at least 60 seconds, not {value}." in flat, flat + + +@pytest.mark.parametrize("command", COMMANDS) +def test_progress_interval_is_in_help(command): + r = _mafutils([command, "--help"]) + assert r.returncode == 0 and "--progress-interval" in r.stdout + + +@pytest.mark.parametrize("command", ["info", "validate"]) +def test_commands_without_progress_do_not_take_the_option(command): + r = _mafutils([command, "--help"]) + assert r.returncode == 0 and "--progress-interval" not in r.stdout + + +@pytest.mark.parametrize("mode_args", [ + ["-b", "id", "-p", "2"], + ["-m", "scaffold", "-i", MAF_FILE + ".scaffold.idx", "-p", "2"], +]) +def test_fetch_outputs_identical_with_progress_off(mode_args, tmp_path): + """Off only removes bookkeeping; the extracted alignments must not change.""" + bed = os.path.join(TEST_DIR, "example.bed") + if "scaffold" in mode_args: + bed = os.path.join(str(tmp_path), "s.bed") + with open(bed, "w") as fp: + fp.write("chr1\t0\t50\tchr1\nchr4\t0\t50\tchr4\n") + trees = {} + for interval in ("60", "0"): + out = os.path.join(str(tmp_path), f"out{interval}") + r = _fetch([MAF_FILE, bed, "-o", out, "--progress-interval", interval] + mode_args) + assert r.returncode == 0, r.stdout + r.stderr + trees[interval] = { + name: open(os.path.join(out, name), "rb").read() + for name in sorted(os.listdir(out)) if name != "maf_fetch.log" + } + assert trees["60"] == trees["0"] and trees["0"] + + +@pytest.mark.parametrize("parallel", [False, True]) +def test_gc_outputs_identical_with_progress_off(parallel, tmp_path): + results = {} + for interval in ("60", "0"): + prefix = os.path.join(str(tmp_path), f"g{interval}") + args = ["gc", MAF_FILE, "-o", prefix, "--progress-interval", interval] + if parallel: + args[2:2] = [MAF_FILE + ".block.idx"] + args += ["-p", "2", "--chunk-size", "3"] + r = _mafutils(args) + assert r.returncode == 0, r.stdout + r.stderr + assert ("PARALLEL" in r.stdout) == parallel + results[interval] = open(prefix + ".gc.csv").read() + assert results["60"] == results["0"] diff --git a/tests/test_scaffold_mode.py b/tests/test_scaffold_mode.py index e77c5b4..57aca0e 100644 --- a/tests/test_scaffold_mode.py +++ b/tests/test_scaffold_mode.py @@ -13,7 +13,7 @@ MAF_FILE = os.path.join(TEST_DIR, "example.maf") SCAFFOLD_INDEX = os.path.join(TEST_DIR, "example.maf.scaffold.idx") -HEADER = "# mafutils-index format=2 maf=t.maf compression=none size=776 mtime=1.0 hash=md5:x\n" +HEADER = "# mafutils-index format=2 maf=t.maf compression=none size=777 mtime=1.0 hash=md5:x\n" def run_fetch(args): diff --git a/tests/test_species_scaffolds.py b/tests/test_species_scaffolds.py new file mode 100644 index 0000000..1a51e49 --- /dev/null +++ b/tests/test_species_scaffolds.py @@ -0,0 +1,273 @@ +import csv +import os +import re +import subprocess +import sys + +import pytest + +from mafutils.stats import SPECIES_SCAFFOLD_HEADERS, writeSpeciesScaffolds + +TEST_DIR = os.path.dirname(__file__) +REPO_ROOT = os.path.abspath(os.path.join(TEST_DIR, "..")) + +MAF_FILE = os.path.join(TEST_DIR, "example.maf") +TABS_MAF = os.path.join(TEST_DIR, "example-tabs.maf") +INDEX_FILE = os.path.join(TEST_DIR, "example.maf.block.idx") + +# Hand-derived from tests/example.maf's s lines: +# (species, scaffold): (src_size, blocks, aligned_bases) +EXPECTED_SCAFFOLDS = { + ("human", "chr1"): (1000, 3, 25), + ("human", "chr2"): (2000, 1, 6), + ("human", "chr3"): (1000, 2, 9), + ("human", "chr4"): (1000, 3, 14), + ("human", "chrX"): (500, 1, 4), + ("chimp", "chr1"): (1000, 3, 20), + ("chimp", "chr2"): (2000, 1, 6), + ("chimp", "chr3"): (1000, 2, 7), + ("chimp", "chr4"): (1000, 2, 8), + ("gorilla", "chr4"): (1000, 1, 4), + ("gorilla", "chrX"): (800, 1, 4), +} +# species -> (scaffolds_in_maf, scaffold_bases_in_maf) +EXPECTED_ROLLUP = {"human": (5, 5500), "chimp": (4, 5000), "gorilla": (2, 1800)} + +ORIGINAL_SPECIES_COLUMNS = [ + "species", "blocks_present", "blocks_missing", "presence_fraction", "copy_lines", + "duplicated_blocks", "duplicate_copies_total", "gaps_total", "nongaps_total", + "gap_fraction_cells", "avg_gaps_per_present_block", + "avg_gap_fraction_per_present_block", "all_gap_blocks", +] + + +def run_stats(args): + return subprocess.run( + [sys.executable, "-m", "mafutils", "stats"] + args, + check=False, cwd=REPO_ROOT, capture_output=True, text=True, + ) + + +def read_tsv(path): + with open(path, encoding="utf-8") as fp: + return list(csv.DictReader(fp, delimiter="\t")) + + +def stats_outputs(tmp_path, maf=MAF_FILE, extra=()): + prefix = os.path.join(str(tmp_path), "out") + args = [maf, "-o", prefix] + list(extra) + if maf == MAF_FILE: + args[1:1] = [INDEX_FILE] + r = run_stats(args) + assert r.returncode == 0, r.stdout + r.stderr + return prefix + + +# --------------------------------------------------------------------------- +# species_scaffolds.tsv +# --------------------------------------------------------------------------- + + +def test_scaffold_table_matches_hand_derived_values(tmp_path): + prefix = stats_outputs(tmp_path) + rows = read_tsv(prefix + ".species_scaffolds.tsv") + got = {(r["species"], r["scaffold"]): (int(r["src_size"]), int(r["blocks"]), int(r["aligned_bases"])) + for r in rows} + assert got == EXPECTED_SCAFFOLDS + + +def test_same_named_scaffolds_in_different_species_are_kept_apart(tmp_path): + """ + human.chr4 and chimp.chr4 are different sequences that share a name; keying + by scaffold alone would merge their blocks and aligned bases. + """ + prefix = stats_outputs(tmp_path) + rows = {(r["species"], r["scaffold"]): r for r in read_tsv(prefix + ".species_scaffolds.tsv")} + assert rows[("human", "chr4")]["aligned_bases"] == "14" + assert rows[("chimp", "chr4")]["aligned_bases"] == "8" + assert rows[("gorilla", "chr4")]["aligned_bases"] == "4" + + +def test_scaffold_table_sorted_by_species_then_largest_first(tmp_path): + prefix = stats_outputs(tmp_path) + rows = read_tsv(prefix + ".species_scaffolds.tsv") + keys = [(r["species"], -int(r["src_size"]), r["scaffold"]) for r in rows] + assert keys == sorted(keys) + human = [r["scaffold"] for r in rows if r["species"] == "human"] + assert human[0] == "chr2" # 2000, the largest + + +def test_aligned_fraction_is_aligned_over_src_size(tmp_path): + prefix = stats_outputs(tmp_path) + for r in read_tsv(prefix + ".species_scaffolds.tsv"): + assert float(r["aligned_fraction"]) == pytest.approx(int(r["aligned_bases"]) / int(r["src_size"])) + + +def test_aligned_fraction_is_not_clamped(tmp_path): + """A region aligned into several blocks can exceed its scaffold length; + that is information, so it must not be capped at 1.0.""" + out = os.path.join(str(tmp_path), "s.tsv") + writeSpeciesScaffolds(out, {("sp", "s1"): {"src_size": 10, "blocks": 3, "aligned_bases": 25}}) + row = read_tsv(out)[0] + assert list(row) == SPECIES_SCAFFOLD_HEADERS + assert row["aligned_fraction"] == "2.50000000" + + +def test_missing_scaffold_name_and_size_are_written_as_dots(tmp_path): + out = os.path.join(str(tmp_path), "s.tsv") + writeSpeciesScaffolds(out, {("sp", ""): {"src_size": None, "blocks": 1, "aligned_bases": 4}}) + row = read_tsv(out)[0] + assert (row["scaffold"], row["src_size"], row["aligned_fraction"]) == (".", ".", ".") + + +def test_tab_delimited_maf_gives_identical_scaffold_table(tmp_path): + os.makedirs(tmp_path / "space") + os.makedirs(tmp_path / "tabs") + space = stats_outputs(tmp_path / "space") + tabs = stats_outputs(tmp_path / "tabs", maf=TABS_MAF) + assert read_tsv(space + ".species_scaffolds.tsv") == read_tsv(tabs + ".species_scaffolds.tsv") + + +# --------------------------------------------------------------------------- +# species.tsv roll-up columns +# --------------------------------------------------------------------------- + + +def test_species_rollup_columns(tmp_path): + prefix = stats_outputs(tmp_path) + rows = {r["species"]: r for r in read_tsv(prefix + ".species.tsv")} + for sp, (n, bases) in EXPECTED_ROLLUP.items(): + assert int(rows[sp]["scaffolds_in_maf"]) == n, sp + assert int(rows[sp]["scaffold_bases_in_maf"]) == bases, sp + + +def test_new_species_columns_are_appended_not_inserted(tmp_path): + """Positional parsers of species.tsv must keep working.""" + prefix = stats_outputs(tmp_path) + header = open(prefix + ".species.tsv").readline().rstrip("\n").split("\t") + assert header[:len(ORIGINAL_SPECIES_COLUMNS)] == ORIGINAL_SPECIES_COLUMNS + assert header[len(ORIGINAL_SPECIES_COLUMNS):] == ["scaffolds_in_maf", "scaffold_bases_in_maf"] + + +def test_expected_but_absent_species_gets_zero_scaffolds(tmp_path): + species_file = os.path.join(str(tmp_path), "sp.txt") + open(species_file, "w").write("human\nchimp\ngorilla\norangutan\n") + prefix = stats_outputs(tmp_path, extra=["--expected-species-file", species_file]) + rows = {r["species"]: r for r in read_tsv(prefix + ".species.tsv")} + assert rows["orangutan"]["scaffolds_in_maf"] == "0" + assert rows["orangutan"]["scaffold_bases_in_maf"] == "0" + + +# --------------------------------------------------------------------------- +# Invariants against existing, independently pinned statistics +# --------------------------------------------------------------------------- + + +def test_per_scaffold_aligned_bases_sum_to_nongaps_total(tmp_path): + """ + Every non-gap base of a species comes from exactly one of its scaffolds, so + the per-scaffold sums must equal the species' existing nongaps_total -- + an independent check on the new accounting. + """ + prefix = stats_outputs(tmp_path) + species = {r["species"]: r for r in read_tsv(prefix + ".species.tsv")} + totals = {} + for r in read_tsv(prefix + ".species_scaffolds.tsv"): + totals[r["species"]] = totals.get(r["species"], 0) + int(r["aligned_bases"]) + for sp, total in totals.items(): + assert total == int(species[sp]["nongaps_total"]), sp + + +def test_per_scaffold_blocks_cover_blocks_present(tmp_path): + """ + sum(blocks) over a species' scaffolds is >= its blocks_present, with + equality when no block carries two copies of the species (true here). + """ + prefix = stats_outputs(tmp_path) + species = {r["species"]: r for r in read_tsv(prefix + ".species.tsv")} + totals = {} + for r in read_tsv(prefix + ".species_scaffolds.tsv"): + totals[r["species"]] = totals.get(r["species"], 0) + int(r["blocks"]) + for sp, total in totals.items(): + assert total == int(species[sp]["blocks_present"]), sp + + +@pytest.mark.parametrize("processes,chunk_size", [(1, 1), (1, 3), (2, 1), (2, 2), (4, 3)]) +def test_scaffold_outputs_invariant_to_chunking(processes, chunk_size, tmp_path): + os.makedirs(tmp_path / "base") + os.makedirs(tmp_path / "run") + base = stats_outputs(tmp_path / "base", extra=["-p", "1", "--chunk-size", "10000"]) + run = stats_outputs(tmp_path / "run", extra=["-p", str(processes), "--chunk-size", str(chunk_size)]) + assert read_tsv(run + ".species_scaffolds.tsv") == read_tsv(base + ".species_scaffolds.tsv") + cols = ["species", "scaffolds_in_maf", "scaffold_bases_in_maf"] + pick = lambda path: [{c: r[c] for c in cols} for r in read_tsv(path)] + assert pick(run + ".species.tsv") == pick(base + ".species.tsv") + + +# --------------------------------------------------------------------------- +# srcSize integrity +# --------------------------------------------------------------------------- + + +def _altered_copy(tmp_path, old, new): + """A copy of example.maf with one s line changed, freshly indexed.""" + text = open(MAF_FILE).read() + assert text.count(old) == 1 + maf = os.path.join(str(tmp_path), "altered.maf") + open(maf, "w").write(text.replace(old, new)) + r = subprocess.run([sys.executable, "-m", "mafutils", "index", maf], + check=False, cwd=REPO_ROOT, capture_output=True, text=True) + assert r.returncode == 0, r.stderr + return maf + + +@pytest.mark.parametrize("processes,chunk_size", [(1, 10000), (2, 3)]) +def test_conflicting_src_size_fails_cleanly(processes, chunk_size, tmp_path): + """ + One species.scaffold reporting two lengths is an inconsistent MAF. It must + fail -- both within a chunk and when the two lines land in different + chunks merged across workers -- naming the scaffold, without a traceback. + """ + maf = _altered_copy(tmp_path, "s chimp.chr4 28 5 + 1000 TGG-C-T", "s chimp.chr4 28 5 + 998 TGG-C-T") + r = run_stats([maf, "-o", os.path.join(str(tmp_path), "out"), + "-p", str(processes), "--chunk-size", str(chunk_size)]) + combined = r.stdout + r.stderr + assert r.returncode != 0 + assert re.search(r"chimp\.chr4 has conflicting srcSize values \((1000 and 998|998 and 1000)\)", combined), combined + assert "Traceback" not in combined + + +def test_unparseable_src_size_affects_only_the_scaffold_table(tmp_path): + maf = _altered_copy(tmp_path, "s chimp.chr4 28 5 + 1000 TGG-C-T", "s chimp.chr4 28 5 + abc TGG-C-T") + prefix = os.path.join(str(tmp_path), "out") + r = run_stats([maf, "-o", prefix]) + assert r.returncode == 0, r.stdout + r.stderr + assert "did not parse as an integer" in (r.stdout + r.stderr) + + overall = {row["metric"]: row["value"] for row in read_tsv(prefix + ".overall.tsv")} + assert overall["scaffold_size_parse_errors"] == "1" + # the line still counts toward existing statistics, unchanged + assert overall["total_sequence_lines"] == "20" + species = {r["species"]: r for r in read_tsv(prefix + ".species.tsv")} + assert species["chimp"]["nongaps_total"] == "41" + # the other chimp.chr4 line (block 8) still supplies its size + rows = {(r["species"], r["scaffold"]): r for r in read_tsv(prefix + ".species_scaffolds.tsv")} + assert rows[("chimp", "chr4")]["src_size"] == "1000" + + +# --------------------------------------------------------------------------- +# Dashboard +# --------------------------------------------------------------------------- + + +def test_dashboard_states_truncation_and_shows_largest_first(tmp_path): + pytest.importorskip("matplotlib") + prefix = stats_outputs(tmp_path, extra=["--html-dashboard", "--dashboard-top-scaffolds", "2"]) + page = open(prefix + ".dashboard.html").read() + summaries = [re.sub(r"<[^>]+>", "", m) for m in re.findall(r".*?", page)] + assert any("human" in s and "showing 2 of 5 scaffolds (largest by length)" in s for s in summaries), summaries + assert any("gorilla" in s and "showing 2 of 2 scaffolds" in s and "largest" not in s for s in summaries) + human = re.search(r"human.*?
    ", page, re.S).group(0) + first = re.search(r"([^<]+)", human).group(1) + assert first == "chr2" + assert "present in the" in page and "MAF" in page From 5d66828db251429c4e33de205acebd567698c814 Mon Sep 17 00:00:00 2001 From: gwct Date: Thu, 22 Oct 2026 10:44:19 -0400 Subject: [PATCH 2/2] Strip ANSI colour from typer output in progress-interval tests --- tests/test_progress.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/tests/test_progress.py b/tests/test_progress.py index 2bb8872..5679687 100644 --- a/tests/test_progress.py +++ b/tests/test_progress.py @@ -2,6 +2,7 @@ import io import logging import os +import re import subprocess import sys @@ -333,6 +334,20 @@ def _mafutils(args): check=False, cwd=REPO_ROOT, capture_output=True, text=True) +ANSI_ESCAPE = re.compile(r"\x1b\[[0-9;]*m") + + +def _plain(text): + """ + typer's own help and error output, as plain text. It is coloured whenever + rich thinks it is on a terminal -- including on GitHub Actions, where + `--progress-interval` arrives as "-", "-progress", "-interval" with escape + codes between them -- and errors are drawn in a box and wrapped. Strip the + codes and the box, and join the lines. + """ + return " ".join(ANSI_ESCAPE.sub("", text).replace("\u2502", " ").split()) + + COMMANDS = ["index", "stats", "gc", "fetch"] @@ -355,22 +370,21 @@ def test_progress_interval_below_minimum_is_a_usage_error(command, value, tmp_pa """Rejected, never clamped: a silently changed value would hide the mistake.""" r = _mafutils([command] + _command_args(command, tmp_path) + ["--progress-interval", value]) assert r.returncode == 2 - # typer draws the error in a box and wraps it; compare with the box and - # line breaks collapsed - flat = " ".join((r.stdout + r.stderr).replace("\u2502", " ").split()) + flat = _plain(r.stdout + r.stderr) assert f"must be 0 (off) or at least 60 seconds, not {value}." in flat, flat @pytest.mark.parametrize("command", COMMANDS) def test_progress_interval_is_in_help(command): r = _mafutils([command, "--help"]) - assert r.returncode == 0 and "--progress-interval" in r.stdout + assert r.returncode == 0 and "--progress-interval" in _plain(r.stdout) @pytest.mark.parametrize("command", ["info", "validate"]) def test_commands_without_progress_do_not_take_the_option(command): r = _mafutils([command, "--help"]) - assert r.returncode == 0 and "--progress-interval" not in r.stdout + # plain text matters here too: in colour, "not in" would pass vacuously + assert r.returncode == 0 and "--progress-interval" not in _plain(r.stdout) @pytest.mark.parametrize("mode_args", [