diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index 2e921ef..a0774bd 100644 --- a/DEVELOPMENT.md +++ b/DEVELOPMENT.md @@ -28,6 +28,30 @@ for user-facing installation/usage docs. missing or headerless), even when the block index matches the MAF file perfectly. +- **The index header carries aggregate counts, and both files of a pair must + get the identical set.** `mafutils index` accumulates blocks / scaffolds / + ref_bases / aln_cols / seq_lines / max_seqs in its existing block loop and + writes them as extra `key=value` tokens, so `mafutils info` can report them + from one line instead of rescanning. Two constraints: + `validate` compares the block and scaffold headers with **exact equality** + (`validate.py`), so block-only keys would make every future validate report + MISMATCH -- the same dict goes into both. And the format stays at **2**: + the change is purely additive, `readIndexHeader` already parses arbitrary + keys, and `tests/test_validate.py` asserts `format == "2"`. + `process_maf_block` returns the integers alongside its stringified row so the + loop never int()s its own output back again; measured cost of the whole + change on a 1.5GB slice, 3 runs each: **4.71s -> 4.77s (~1.3%)**, and a full + 42.7GB rebuild produced aggregates matching an independent awk pass exactly + (`blocks=8094203 scaffolds=61 ref_bases=2728222451 aln_cols=3031569326 + seq_lines=103950448 max_seqs=15`) with all 8.09M rows byte-identical. +- **`mafutils info` never scans the MAF for counts.** Header first; failing + that, stream the *index* (~52 MB/s) with a warning naming its size; failing + that, report only what needs no index. Species are the exception and are + always labelled a sample -- they are not in the index, and collecting them at + 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. + ## Compression / index internals - **`mafutils/lib/bgzf.py` is vendored from Biopython, not a dependency.** diff --git a/README.md b/README.md index 5f60237..b3e1d6d 100644 --- a/README.md +++ b/README.md @@ -3,9 +3,10 @@ `mafutils` is a command-line toolkit for indexing, extracting, and summarizing MAF (Multiple Alignment Format) files. -It currently provides five commands: +It currently provides six commands: - `mafutils index` +- `mafutils info` - `mafutils fetch` - `mafutils stats` - `mafutils gc` @@ -51,6 +52,12 @@ Create block and scaffold indexes for a MAF (defaults to mafutils index input.maf ``` +Print a quick overview of an indexed MAF (sub-second on any size): + +```bash +mafutils info input.maf +``` + Fetch trimmed MAF regions from a BED file (index defaults to `input.maf.block.idx`): @@ -154,6 +161,63 @@ Arguments: | `BLOCK_INDEX` | Output block index path (default: `.block.idx`) | | `SCAFFOLD_INDEX` | Output scaffold index path (default: `.scaffold.idx`) | +### `mafutils info` + +Print a quick overview of a MAF without reading it. + +```bash +mafutils info [OPTIONS] MAF_FILE +``` + +``` +MAF : cricetid-15spec.Mus_musculus.nodupes.maf +size : 42.70 GB (none) +index : matches MAF (size + mtime) + +ref scaffolds : 61 +blocks : 8,094,203 +ref bases : 2,728,222,451 +aln columns : 3,031,569,326 +seq lines : 103,950,448 (all species, including the reference) +max seqs/block : 15 + +species : 15 (sampled from the first 1,000 block(s) -- possibly incomplete) + Arvicola_amphibius, Cricetulus_griseus, Cricetus_cricetus, ... + +For exact per-species presence across the whole file, run `mafutils stats`. +``` + +Measured on that 42.70 GB file: **1.3 s**, or **0.7 s** with `--sample-blocks 0`. + +| Option | Description | +|---|---| +| `--index`, `-i` | Block index (default: `.block.idx`) | +| `--scaffold-index` | Scaffold index (default: `.scaffold.idx`) | +| `--sample-blocks` | Blocks to sample for species names (default 1000); `0` disables the sample and reads no MAF data at all | +| `--all-species` | List every sampled species instead of the first 20 | + +**Where the numbers come from.** `mafutils index` records the block counts in +the index header, so `info` normally reads nothing but that one line. Three +tiers, and the MAF is never scanned for these counts: + +1. Index built by mafutils 0.7 or later — instant, header only. +2. Older index — falls back to streaming the *index* (~52 MB/s, so ~9 s for a + 445 MB index but ~12 min for a 35 GB one), with a warning naming the index + size and recommending a rebuild. +3. No index — reports what needs none and points at `mafutils index`. + +**`species` is a sample, not a count.** Species names are not recorded in the +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. + +**Column meanings.** `ref scaffolds` and `ref bases` are measured on the +*reference* sequence. `seq lines` counts every `s` line including the +reference. `aln columns` is the block width, gap columns included, which all +species share. + ### `mafutils fetch` Fetch regions or scaffolds from a MAF using an existing index. diff --git a/mafutils/cli.py b/mafutils/cli.py index 6a9329e..346c422 100644 --- a/mafutils/cli.py +++ b/mafutils/cli.py @@ -6,6 +6,7 @@ from . import fetch as fetch_mod from . import gc as gc_mod from . import index as index_mod +from . import info as info_mod from . import stats as stats_mod from . import validate as validate_mod @@ -19,6 +20,7 @@ app.command("fetch", help="Fetch alignment blocks from a MAF using a BED file and an existing index.")(fetch_mod.fetch_command) app.command("gc", help="Calculate per-species GC content from a MAF file.")(gc_mod.gc_command) app.command("index", help="Create block and scaffold indexes for a MAF file.")(index_mod.index_command) +app.command("info", help="Print a quick overview of a MAF file from its index.")(info_mod.info_command) app.command("stats", help="Summarize an indexed MAF at overall, species, and block levels.")(stats_mod.stats_command) app.command("validate", help="Check whether a MAF file's index is still trustworthy.")(validate_mod.validate_command) diff --git a/mafutils/index.py b/mafutils/index.py index 871064e..7963b00 100644 --- a/mafutils/index.py +++ b/mafutils/index.py @@ -17,12 +17,22 @@ def process_maf_block(block): + """ + Returns (row, ref_len, aln_cols, num_seqs) for one block: the index row as + strings, plus the three integers the header aggregates need. + + The integers are returned rather than re-parsed from `row` on purpose -- + this loop runs once per block (233M times on a whole-genome MAF) and is the + hot path the v0.6.0 indexing speedup optimised, so int()-ing its own output + back again would be a needless per-block cost. + """ ref_seq = block[1].split() ref_scaff = ref_seq[1].split(".", 1)[1] line_len = str(len(block[1])) - num_seqs = str(len(block) - 1) - seq_len = str(len(ref_seq[6])) - return [ref_scaff, ref_seq[2], ref_seq[3], seq_len, line_len, num_seqs] + num_seqs = len(block) - 1 # every s line, including the reference + seq_len = len(ref_seq[6]) # block width, gap columns included + row = [ref_scaff, ref_seq[2], ref_seq[3], str(seq_len), line_len, str(num_seqs)] + return row, int(ref_seq[3]), seq_len, num_seqs def run_index(maf_file, block_index_path, scaffold_index_path): @@ -52,10 +62,27 @@ def run_index(maf_file, block_index_path, scaffold_index_path): # (see iterMafBlocks/openMafHashing). use_binary = maf_compression in ("none", "gz") + # Header aggregates, accumulated as we go so `mafutils info` can + # report them later from the header alone rather than rescanning. + n_blocks = 0 + total_ref_bases = 0 + total_aln_cols = 0 + total_seq_lines = 0 + max_seqs = 0 + distinct_scaffolds = set() + for block, block_start, block_end in COMMON.iterMafBlocks(maf_stream, binary=use_binary): - block_info = process_maf_block(block) + block_info, ref_len, aln_cols, num_seqs = process_maf_block(block) ref_scaffold = block_info[0] + n_blocks += 1 + total_ref_bases += ref_len + total_aln_cols += aln_cols + total_seq_lines += num_seqs + if num_seqs > max_seqs: + max_seqs = num_seqs + distinct_scaffolds.add(ref_scaffold) + mdx_line = block_info + [str(block_start), str(block_end)] block_stream.write("\t".join(mdx_line) + "\n") @@ -76,13 +103,24 @@ def run_index(maf_file, block_index_path, scaffold_index_path): content_hash = f"md5:{hash_obj.hexdigest()}" + # Identical in both headers on purpose -- validate cross-checks the two + # for exact equality, and these describe the MAF, not the index file. + aggregates = { + "blocks": n_blocks, + "scaffolds": len(distinct_scaffolds), + "ref_bases": total_ref_bases, + "aln_cols": total_aln_cols, + "seq_lines": total_seq_lines, + "max_seqs": max_seqs, + } + with open(block_index_path, "w") as block_stream: - COMMON.writeIndexHeader(block_stream, maf_file, maf_compression, size, mtime, content_hash) + COMMON.writeIndexHeader(block_stream, maf_file, maf_compression, size, mtime, content_hash, aggregates) with open(block_tmp_path, "r") as tmp_fp: shutil.copyfileobj(tmp_fp, block_stream) with open(scaffold_index_path, "w") as scaffold_stream: - COMMON.writeIndexHeader(scaffold_stream, maf_file, maf_compression, size, mtime, content_hash) + COMMON.writeIndexHeader(scaffold_stream, maf_file, maf_compression, size, mtime, content_hash, aggregates) with open(scaffold_tmp_path, "r") as tmp_fp: shutil.copyfileobj(tmp_fp, scaffold_stream) diff --git a/mafutils/info.py b/mafutils/info.py new file mode 100644 index 0000000..8b2bc2d --- /dev/null +++ b/mafutils/info.py @@ -0,0 +1,270 @@ +############################################################################# +# Print a quick overview of a MAF file from its index. +# +# Gregg Thomas + Claude, September 2026 +############################################################################# + +""" +mafutils info + +Prints a fast overview of a MAF: size, compression, index freshness, and the +block/scaffold counts recorded in the index header. + +Deliberately cheap. Three tiers, in order of cost, and the MAF itself is never +scanned for these numbers: + + 1. Index header carries the aggregate keys (built by mafutils >= 0.7): + header read only, instant regardless of MAF size. + 2. Keys absent (older index): fall back to streaming the *index*, which is + orders of magnitude smaller than the MAF -- ~52 MB/s, so ~9s for a 445MB + index but ~12 min for a 35GB one. Warns up front with the index size and + recommends rebuilding, so a long wait is never a surprise. + 3. No index: report only what needs no index and point at `mafutils index`. + +Species are SAMPLED from the first N blocks (default 1000) and always labelled +as such: species names are not recorded in the index, and collecting them +during indexing would parse every non-reference `s` line, measured at 2.06x +slower indexing -- undoing most of the speedup in v0.6.0. For an exhaustive +per-species picture use `mafutils stats`, which writes .species.tsv. +""" + +import logging +import os +import sys +from typing import Annotated, Optional + +import typer + +from mafutils.lib import common as COMMON + +DEFAULT_SAMPLE_BLOCKS = 1000 +SPECIES_PREVIEW = 20 + +############################################################################# + + +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 + index predates them. All-or-nothing on purpose: a partial set would mean + reporting some counts and silently omitting others. + """ + if not header: + return None + if not all(key in header for key in COMMON.INDEX_AGGREGATE_KEYS): + return None + try: + return {key: int(header[key]) for key in COMMON.INDEX_AGGREGATE_KEYS} + except ValueError: + return None + + +def scanIndexAggregates(index_file): + """ + Recomputes the aggregates by streaming the block index. Only used for + indexes built before the header recorded them. + + Streams and accumulates -- never materializes the rows. A whole-genome + index is 35GB and ~233M rows, so holding them would repeat the OOM this + project already fixed in gc/stats. + """ + blocks = ref_bases = aln_cols = seq_lines = max_seqs = 0 + scaffolds = set() + with open(index_file, "r", encoding="utf-8") as fp: + for line in fp: + if line.startswith("#") or not line.strip(): + continue + fields = line.rstrip("\n").split("\t") + if len(fields) < 8: + continue + try: + ref_len = int(fields[2]) + cols = int(fields[3]) + nseq = int(fields[5]) + except ValueError: + continue + blocks += 1 + ref_bases += ref_len + aln_cols += cols + seq_lines += nseq + if nseq > max_seqs: + max_seqs = nseq + scaffolds.add(fields[0]) + return { + "blocks": blocks, + "scaffolds": len(scaffolds), + "ref_bases": ref_bases, + "aln_cols": aln_cols, + "seq_lines": seq_lines, + "max_seqs": max_seqs, + } + + +def sampleSpecies(maf_file, compression, max_blocks): + """ + Collects species names from the first `max_blocks` blocks. Returns + (sorted names, blocks actually seen). + + Sequence lines are detected with split()[0] == "s", never + startswith("s ") -- real MAFs are tab-delimited and the space form matches + nothing on them, a false negative this project has already been bitten by. + """ + species = set() + blocks = 0 + with COMMON.openMaf(maf_file, compression, "rt") as fp: + for line in fp: + if not line.strip(): + continue + fields = line.split() + if fields[0] == "a": + blocks += 1 + if blocks > max_blocks: + blocks -= 1 + break + elif fields[0] == "s" and len(fields) > 1: + # src is "species.scaffold"; the species is everything before + # the first dot (mirrors fetch.speciesFromSrc, which lives in + # fetch.py rather than common). + species.add(fields[1].split(".", 1)[0]) + return sorted(species), blocks + + +############################################################################# + + +def run_info(args, LOG): + maf_compression = COMMON.detectCompression(args.maf_file) + if maf_compression not in ("none", "gz", "bgzip"): + LOG.error(f"Unsupported MAF compression: {maf_compression}") + sys.exit(1) + + index_file = args.index_file or COMMON.deriveBlockIndexPath(args.maf_file) + scaffold_index_file = args.scaffold_index_file or COMMON.deriveScaffoldIndexPath(args.maf_file) + have_index = os.path.isfile(index_file) + + header = COMMON.readIndexHeader(index_file) if have_index else None + aggregates = readAggregates(header) + + if have_index and aggregates is None: + index_size = os.path.getsize(index_file) + 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), + ) + aggregates = scanIndexAggregates(index_file) + + # Scaffold names come from the scaffold index, which is tiny (a couple of + # KB even for a whole genome) -- worth reading whenever it exists. + scaffold_names = [] + if os.path.isfile(scaffold_index_file): + seen = [] + with open(scaffold_index_file, "r", encoding="utf-8") as fp: + for line in fp: + if line.startswith("#") or not line.strip(): + continue + name = line.split("\t", 1)[0] + if name not in seen: + seen.append(name) + scaffold_names = seen + + species, sampled_blocks = ([], 0) + if args.sample_blocks > 0: + species, sampled_blocks = sampleSpecies(args.maf_file, maf_compression, args.sample_blocks) + + # ---- output ------------------------------------------------------- + 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})") + + if not have_index: + out.append(f"{'index':<15}: not found at {index_file}") + else: + comparison = COMMON.compareIndexHeader(header, args.maf_file, maf_compression) + problems = [k for k, v in comparison.items() if v is not None and not v["match"]] + if header is None: + state = "present, but has no mafutils header (built by an older version)" + elif problems: + state = "STALE -- " + ", ".join(f"{k} differs" for k in problems) + else: + state = "matches MAF (size + mtime)" + out.append(f"{'index':<15}: {state}") + + out.append("") + if scaffold_names: + out.append(f"{'ref scaffolds':<15}: {len(scaffold_names):,}") + elif aggregates: + out.append(f"{'ref scaffolds':<15}: {aggregates['scaffolds']:,}") + + if aggregates: + out.append(f"{'blocks':<15}: {aggregates['blocks']:,}") + out.append(f"{'ref bases':<15}: {aggregates['ref_bases']:,}") + out.append(f"{'aln columns':<15}: {aggregates['aln_cols']:,}") + out.append(f"{'seq lines':<15}: {aggregates['seq_lines']:,} (all species, including the reference)") + out.append(f"{'max seqs/block':<15}: {aggregates['max_seqs']:,}") + else: + out.append(f"{'blocks':<15}: unavailable without an index -- run `mafutils index` first") + + if args.sample_blocks > 0: + out.append("") + out.append( + f"{'species':<15}: {len(species):,} " + f"(sampled from the first {sampled_blocks:,} block(s) -- possibly incomplete)" + ) + shown = species[:SPECIES_PREVIEW] + line = " " + for name in shown: + if len(line) + len(name) + 2 > 78: + out.append(line.rstrip()) + line = " " + line += name + ", " + out.append(line.rstrip().rstrip(",")) + if len(species) > SPECIES_PREVIEW: + out.append(f" ... and {len(species) - SPECIES_PREVIEW:,} more (--all-species to list them)") + out.append("") + out.append("For exact per-species presence across the whole file, run `mafutils stats`.") + + print("\n".join(out)) + + +############################################################################# + + +def info_command( + maf_file: Annotated[str, typer.Argument(help="Input MAF file (.maf, .maf.gz, or bgzip-compressed .maf)")], + index_file: Annotated[Optional[str], typer.Option("--index", "-i", help="Block index (default: .block.idx)")] = None, + scaffold_index_file: Annotated[Optional[str], typer.Option("--scaffold-index", help="Scaffold index (default: .scaffold.idx)")] = None, + sample_blocks: Annotated[int, typer.Option("--sample-blocks", help="Blocks to sample for species names; 0 disables the sample and reads no MAF data at all.")] = DEFAULT_SAMPLE_BLOCKS, + all_species: Annotated[bool, typer.Option("--all-species", help="List every sampled species name instead of the first 20.")] = False, +) -> None: + if sample_blocks < 0: + print("--sample-blocks must be >= 0", file=sys.stderr) + raise typer.Exit(1) + + global SPECIES_PREVIEW + if all_species: + SPECIES_PREVIEW = sys.maxsize + + logging.basicConfig(format="[ %(levelname)s ] %(message)s", level=logging.INFO) + LOG = logging.getLogger("maf_info_logger") + + from types import SimpleNamespace + args = SimpleNamespace( + maf_file=maf_file, + index_file=index_file, + scaffold_index_file=scaffold_index_file, + sample_blocks=sample_blocks, + ) + run_info(args, LOG) diff --git a/mafutils/lib/common.py b/mafutils/lib/common.py index 57fdabc..8d86cf5 100644 --- a/mafutils/lib/common.py +++ b/mafutils/lib/common.py @@ -287,10 +287,37 @@ def deriveScaffoldIndexPath(maf_file): INDEX_HEADER_PREFIX = "# mafutils-index" -def writeIndexHeader(stream, maf_file, compression, size, mtime, content_hash): +# Aggregate counts recorded in the index header at build time, so `mafutils +# info` can report them without reading anything but the header. Computed for +# free while indexing (it already streams every block); an index that predates +# them simply lacks the keys and consumers fall back to scanning. +INDEX_AGGREGATE_KEYS = ("blocks", "scaffolds", "ref_bases", "aln_cols", "seq_lines", "max_seqs") + + +def writeIndexHeader(stream, maf_file, compression, size, mtime, content_hash, aggregates=None): + """ + Writes the one-line index header. + + `aggregates` is an optional dict of INDEX_AGGREGATE_KEYS. The SAME values + must be written to both the block and scaffold index of a pair: `mafutils + validate` cross-checks the two headers with exact equality, so keys present + in only one would make every subsequent validate report MISMATCH. They + describe the MAF rather than the index file, so identical values in both is + also the semantically right thing. + + Kept at format=2 deliberately: extra key=value tokens are additive and + readIndexHeader already parses arbitrary keys, with every consumer treating + a missing key as "can't check this dimension". Bumping the format would + imply an incompatibility that does not exist. + """ + extra = "" + if aggregates: + extra = "".join( + f" {key}={aggregates[key]}" for key in INDEX_AGGREGATE_KEYS if key in aggregates + ) stream.write( f"{INDEX_HEADER_PREFIX} format=2 maf={os.path.basename(maf_file)} " - f"compression={compression} size={size} mtime={mtime} hash={content_hash}\n" + f"compression={compression} size={size} mtime={mtime} hash={content_hash}{extra}\n" ) diff --git a/tests/test_info.py b/tests/test_info.py new file mode 100644 index 0000000..00544f9 --- /dev/null +++ b/tests/test_info.py @@ -0,0 +1,196 @@ +import os +import re +import subprocess +import sys + +import pytest + +from mafutils.lib import common as COMMON + +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") + +# Hand-verifiable against tests/example.maf, and cross-checked against the +# independent values already pinned in test_stats.py (10 blocks, 61 alignment +# columns, 20 sequence lines). +EXPECTED = { + "blocks": 10, + "scaffolds": 5, # chr1..chr4, chrX -- chr4 appears in two runs, one scaffold + "ref_bases": 59, + "aln_cols": 61, + "seq_lines": 20, + "max_seqs": 2, +} + + +def run_info(args): + return subprocess.run( + [sys.executable, "-m", "mafutils", "info"] + args, + check=False, cwd=REPO_ROOT, capture_output=True, text=True, + ) + + +def build_index(tmp_path, maf=MAF_FILE, name="t.maf"): + """A freshly built index, i.e. one whose header carries the aggregates.""" + import shutil + copy = os.path.join(str(tmp_path), name) + shutil.copyfile(maf, copy) + r = subprocess.run([sys.executable, "-m", "mafutils", "index", copy], + check=False, cwd=REPO_ROOT, capture_output=True, text=True) + assert r.returncode == 0, r.stderr + return copy + + +def value(output, label): + m = re.search(rf"^{re.escape(label)}\s*:\s*([0-9,]+)", output, re.M) + assert m, f"{label!r} not found in:\n{output}" + return int(m.group(1).replace(",", "")) + + +# --------------------------------------------------------------------------- +# Index header aggregates +# --------------------------------------------------------------------------- + + +def test_index_records_aggregates_in_the_header(tmp_path): + maf = build_index(tmp_path) + header = COMMON.readIndexHeader(COMMON.deriveBlockIndexPath(maf)) + for key, expected in EXPECTED.items(): + assert int(header[key]) == expected, key + + +def test_both_index_headers_are_identical(tmp_path): + """ + validate cross-checks the two headers with exact equality, so aggregates + present in only one would make every future validate report MISMATCH. + """ + maf = build_index(tmp_path) + block = open(COMMON.deriveBlockIndexPath(maf)).readline() + scaffold = open(COMMON.deriveScaffoldIndexPath(maf)).readline() + assert block == scaffold + + +def test_validate_still_verifies_with_aggregates(tmp_path): + maf = build_index(tmp_path) + r = subprocess.run([sys.executable, "-m", "mafutils", "validate", maf], + check=False, cwd=REPO_ROOT, capture_output=True, text=True) + assert r.returncode == 0, r.stdout + r.stderr + assert "VERIFIED" in (r.stdout + r.stderr) + + +def test_format_version_unchanged(tmp_path): + """Additive keys only -- bumping the format would imply an incompatibility + that does not exist, and old readers already ignore unknown keys.""" + maf = build_index(tmp_path) + header = COMMON.readIndexHeader(COMMON.deriveBlockIndexPath(maf)) + assert header["format"] == "2" + + +# --------------------------------------------------------------------------- +# info: the three tiers +# --------------------------------------------------------------------------- + + +def test_info_reports_exact_values_from_the_header(tmp_path): + maf = build_index(tmp_path) + r = run_info([maf]) + assert r.returncode == 0, r.stderr + out = r.stdout + assert value(out, "blocks") == EXPECTED["blocks"] + assert value(out, "ref scaffolds") == EXPECTED["scaffolds"] + assert value(out, "ref bases") == EXPECTED["ref_bases"] + assert value(out, "aln columns") == EXPECTED["aln_cols"] + assert value(out, "seq lines") == EXPECTED["seq_lines"] + assert value(out, "max seqs/block") == EXPECTED["max_seqs"] + # no fallback warning when the header already has them + assert "predates" not in (r.stdout + r.stderr) + + +def test_info_falls_back_to_scanning_an_older_index(tmp_path): + """ + An index without the aggregate keys must still produce the right numbers, + warn that it is scanning, and name the index size so a long wait on a large + index is never a surprise. + """ + maf = build_index(tmp_path) + index_path = COMMON.deriveBlockIndexPath(maf) + lines = open(index_path).readlines() + # strip the aggregate keys, leaving a v0.6.0-shaped header + lines[0] = re.sub(r"\s+(blocks|scaffolds|ref_bases|aln_cols|seq_lines|max_seqs)=\d+", "", lines[0]) + open(index_path, "w").writelines(lines) + assert "blocks=" not in open(index_path).readline() + + r = run_info([maf]) + assert r.returncode == 0, r.stderr + combined = r.stdout + r.stderr + assert "predates the aggregate header fields" in combined + assert "Rebuild with `mafutils index`" in combined + assert re.search(r"scanning the index \([\d.]+ \w+\)", combined), combined + # and the scanned numbers equal what the header would have said + for label, key in (("blocks", "blocks"), ("ref bases", "ref_bases"), + ("aln columns", "aln_cols"), ("seq lines", "seq_lines")): + assert value(r.stdout, label) == EXPECTED[key], label + + +def test_info_without_an_index_says_so_and_does_not_scan(tmp_path): + import shutil + copy = os.path.join(str(tmp_path), "lonely.maf") + shutil.copyfile(MAF_FILE, copy) + + r = run_info([copy]) + assert r.returncode == 0, r.stderr + out = r.stdout + assert "index" in out and "not found" in out + assert "unavailable without an index" in out + assert "run `mafutils index` first" in out + # species still works -- it reads the MAF directly, not the index + assert value(out, "species") == 3 + + +# --------------------------------------------------------------------------- +# Species sampling +# --------------------------------------------------------------------------- + + +def test_species_are_labelled_as_a_sample(tmp_path): + maf = build_index(tmp_path) + out = run_info([maf]).stdout + assert value(out, "species") == 3 + assert "possibly incomplete" in out + assert "sampled from the first" in out + for name in ("human", "chimp", "gorilla"): + assert name in out + assert "mafutils stats" in out + + +def test_sample_blocks_zero_skips_the_species_section(tmp_path): + maf = build_index(tmp_path) + out = run_info([maf, "--sample-blocks", "0"]).stdout + # "species" still appears in the seq-lines note, so assert on the section + assert "possibly incomplete" not in out + assert not re.search(r"^species\s*:", out, re.M), out + assert "mafutils stats" not in out + assert value(out, "blocks") == EXPECTED["blocks"] + + +def test_tab_delimited_maf_yields_the_same_species(tmp_path): + """ + Real MAFs are tab-delimited; detecting sequence lines with + startswith("s ") silently matches nothing on them, a false negative this + project has hit before. + """ + os.makedirs(tmp_path / "tabs", exist_ok=True) + tabs = build_index(tmp_path / "tabs", maf=TABS_MAF, name="tabs.maf") + out = run_info([tabs]).stdout + assert value(out, "species") == 3 + for name in ("human", "chimp", "gorilla"): + assert name in out + + +def test_negative_sample_blocks_is_rejected(tmp_path): + maf = build_index(tmp_path) + r = run_info([maf, "--sample-blocks", "-1"]) + assert r.returncode != 0