From 4af40e9f9796f90ed664300f35718446d1835f74 Mon Sep 17 00:00:00 2001 From: gwct Date: Sun, 8 Nov 2026 21:08:39 -0500 Subject: [PATCH 1/3] Keep fetch's per-element locus summary from growing with region size In MAF mode, fetchByRegion kept one chunk dict per `s` line of a region just to compute the per-element summary columns at the end. A 69 Mb region of a 45-way bird MAF was OOM-killed at 6 GB, where v0.5.0 used 307 MB. ElementLociAccumulator now folds in each chunk as it arrives. A species keeps packed start/size/block arrays only while it is single-scaffold and single-strand. The classification and gap computation are unchanged, and summarizeElementLoci is now a wrapper over the accumulator. FASTA mode holds every chunk anyway, so it still summarizes once at the end, as before. - Worker peak RSS in MAF mode: 93 Mb region 1,531 -> 94 MB, 30 Mb 420 -> 56 MB, 10 Mb 113 -> 44 MB. - Outputs, summary included, are byte-identical to v0.8.1 on a 979,442-region BED in MAF and FASTA mode. Runtime is within noise. - tests/test_loci_accumulator.py checks it against a verbatim copy of v0.8.1's code on 12,000 random elements. Co-Authored-By: Claude Opus 5.5 --- DEVELOPMENT.md | 47 ++++++ mafutils/fetch.py | 200 +++++++++++++++++++------ tests/test_loci_accumulator.py | 260 +++++++++++++++++++++++++++++++++ 3 files changed, 466 insertions(+), 41 deletions(-) create mode 100644 tests/test_loci_accumulator.py diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index 192dfb1..5cf045e 100644 --- a/DEVELOPMENT.md +++ b/DEVELOPMENT.md @@ -155,6 +155,53 @@ for user-facing installation/usage docs. 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. +- **The per-element locus summary no longer grows with region size** + (`ElementLociAccumulator`). + - **The problem:** from v0.6.0, `fetchByRegion` kept one chunk dict per `s` + line of the region (`element_chunks`) just to compute the 12 summary + columns at the end. That is ~450 B per line, so a 69 Mb region of a + 45-way bird MAF (40.5M `s` lines) was OOM-killed at 6 GB, where v0.5.0 + used 307 MB. + - **The fix (MAF output):** each chunk is folded in as it arrives. A + species keeps packed arrays (start/size/block: 16 B per contributing + chunk) only while it is single-scaffold and single-strand, and drops them + as soon as it is multi-scaffold or multi-strand, because those classes use + no gaps. The final computation is the old one unchanged: class + precedence, stable sort by start, boundary attribution and tie-breaking. + `summarizeElementLoci` is now a wrapper over it. + - **Arrays only from the second chunk:** a species' first chunk is kept in + plain fields. Creating three arrays per species up front made the summary + 8% slower than v0.8.1 on small elements (22.5 vs 20.8 µs for a 15-species + element). Lazily it is 18% faster (17.6 vs 21.4 µs), because most species + in a small element have one chunk. + - **FASTA output is unchanged:** it keeps every chunk dict anyway (headers, + `--loci-table`), so it still collects references and calls + `summarizeElementLoci` once at the end, exactly as before. + - **Equivalence:** `tests/test_loci_accumulator.py` checks it against a + verbatim copy of v0.8.1's code on 12,000 random elements (ties, + out-of-order and overlapping chunks, zero sizes, multi + scaffold/strand, repeated or missing `block_index`). + - **Measured on the hamster slice, MAF mode,** worker peak RSS: 93 Mb region + 1,531 -> 94 MB (and 64.6 -> 56.1s); 30 Mb 420 -> 56 MB; 10 Mb 113 -> 44 + MB. Outputs, summary included, are byte-identical to v0.8.1. FASTA mode + must hold every species' sequence for the region, so a huge region stays + large: 93 Mb was 6,719 MB for v0.8.1 and 6,735 MB for the fix. + - **Measuring pitfall:** forkserver workers (Python 3.14) are invisible to + `/usr/bin/time`. Use Python 3.12 (fork), or sample `/proc//status` + VmHWM across the process tree. Also compare both versions run from copies + outside the project directory. Run from the project directory, the same + code peaked ~350 MB (5%) higher on the 93 Mb FASTA region: pymalloc + fragmentation at a 6.7 GB peak shifts with trivial import-time + differences. With `PYTHONMALLOC=malloc` both versions measured + 7,046-7,050 MB. + - **Full hamster BED** (979,442 regions, `-p 8`), v0.8.1 vs fix: all + 979,444 output files byte-identical in MAF and FASTA mode. + - MAF mode: 94.7s vs 94.2s. + - FASTA mode: 119.6s vs 120.4s over six pairs, three in each run order. + Per-pair differences ran -3.8s to +4.7s against a node load of 9-20, + i.e. noise. + - **Not changed:** `block_stats`, one small list per block, about 44 MB at + 93 Mb. - **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`, diff --git a/mafutils/fetch.py b/mafutils/fetch.py index 086870e..5784d2a 100644 --- a/mafutils/fetch.py +++ b/mafutils/fetch.py @@ -77,6 +77,7 @@ from types import SimpleNamespace from typing import Annotated, Optional from concurrent.futures import FIRST_COMPLETED, ProcessPoolExecutor, wait +from array import array from collections import defaultdict, OrderedDict import typer @@ -992,10 +993,11 @@ def chunkGaps(chunks): ] -def summarizeElementLoci(per_species_chunks): +class ElementLociAccumulator: """ Per-ELEMENT summary of how each species' contribution looks, as the values - for ELEMENT_LOCI_HEADERS. + for ELEMENT_LOCI_HEADERS -- built up one chunk at a time with add(). + Classification is irreducibly per-species -- one real element was contiguous for 5 species, split for 6 and empty for 3 -- so this reports @@ -1022,47 +1024,156 @@ def summarizeElementLoci(per_species_chunks): cases) has a gap spanning more than one boundary; it is attributed to the first boundary it crosses, and the per-chunk loci table carries block_index/gap_to_next so that case stays inspectable. + + Memory. This used to keep every chunk of the region as a dict until the + region ended, ~450 B per s line: linear in region size, so a 69 Mb region of + a 45-way MAF (40.5M s lines) was OOM-killed at 6 GB. Now a species keeps + only packed arrays (16 B per contributing chunk, created at its second + chunk), and drops them the moment it becomes multi-scaffold or + multi-strand, because those classes use no gaps (chunkGaps returns []). The values are exactly what the dict version + computed -- same class precedence, same stable sort by start, same boundary + attribution and tie-breaking -- which tests/test_loci_accumulator.py checks + against a verbatim copy of it. """ - counts = {"single": 0, "contiguous": 0, "split": 0, - "multi_scaffold": 0, "multi_strand": 0, "no_bases": 0} - all_gaps = [] - per_species_total = [] - # boundary -> {species: gap}, keyed by the block ordinal the gap starts at - by_boundary = defaultdict(dict) - for species, chunks in per_species_chunks.items(): - locus_class, _max_gap = classifyChunks(chunks) - if locus_class in counts: - counts[locus_class] += 1 + # A chunk with no block_index: its boundary falls back to i + 1, as before. + NO_BLOCK_INDEX = -1 + + def __init__(self): + # species -> [n_chunks, n_contributing, scaffold, strand, + # multi_scaffold, multi_strand, starts, sizes, blocks, + # first_start, first_size, first_block] + # in first-appearance order. A species' first contributing chunk is + # kept in the last three fields; the arrays are created only when a + # second one arrives (most species in a small element have one chunk, + # so they never allocate any), and are dropped again once the species + # is multi-scaffold or multi-strand. + self.species = {} + + def addSpecies(self, species): + rec = self.species.get(species) + if rec is None: + rec = self.species[species] = [0, 0, None, None, False, False, None, None, None, 0, 0, 0] + return rec + + def add(self, species, scaffold, strand, start, size, block_index): + """One chunk. Zero-size chunks only register the species (n.species).""" + rec = self.species.get(species) + if rec is None: + rec = self.addSpecies(species) + rec[0] += 1 + if size <= 0: + return + if block_index is None: + block_index = self.NO_BLOCK_INDEX + n_contrib = rec[1] + rec[1] = n_contrib + 1 + if n_contrib == 0: + rec[2] = scaffold + rec[3] = strand + rec[9] = start + rec[10] = size + rec[11] = block_index + return + if rec[4]: + return # multi_scaffold: final + if scaffold != rec[2]: + rec[4] = True # outranks multi_strand + rec[6] = rec[7] = rec[8] = None + return + if rec[5]: + return # multi_strand, one scaffold so far + if strand != rec[3]: + rec[5] = True + rec[6] = rec[7] = rec[8] = None + return + starts = rec[6] + if starts is None: # second chunk: start the arrays + starts = rec[6] = array('q', (rec[9],)) + rec[7] = array('i', (rec[10],)) + rec[8] = array('i', (rec[11],)) + starts.append(start) + rec[7].append(size) + rec[8].append(block_index) + + def summary(self): + counts = {"single": 0, "contiguous": 0, "split": 0, + "multi_scaffold": 0, "multi_strand": 0, "no_bases": 0} + all_gaps = [] + per_species_total = [] + # boundary -> {species: gap}, keyed by the block ordinal the gap starts at + by_boundary = defaultdict(dict) + + for species, rec in self.species.items(): + n_chunks, n_contrib, _scaffold, _strand, multi_scaffold, multi_strand, starts, sizes, blocks = rec[:9] + if n_chunks == 0: + continue # registered with no chunks: no class + if n_contrib == 0: + counts["no_bases"] += 1 + continue + if n_contrib == 1: + counts["single"] += 1 + continue + if multi_scaffold: + counts["multi_scaffold"] += 1 + continue + if multi_strand: + counts["multi_strand"] += 1 + continue - gaps = chunkGaps(chunks) - nonzero = [g for g in gaps if g != 0] - if nonzero: + # sorted() is stable, so equal starts keep insertion order, as + # sorting the chunk dicts by start did. + order = sorted(range(n_contrib), key=starts.__getitem__) + gaps = [ + starts[order[i + 1]] - (starts[order[i]] + sizes[order[i]]) + for i in range(n_contrib - 1) + ] + nonzero = [g for g in gaps if g != 0] + if not nonzero: + counts["contiguous"] += 1 + continue + counts["split"] += 1 all_gaps.extend(nonzero) per_species_total.append(sum(nonzero)) - contributing = sorted((c for c in chunks if c['size'] > 0), key=lambda c: c['start']) for i, gap in enumerate(gaps): if gap != 0: - by_boundary[contributing[i].get('block_index', i + 1)][species] = gap + block = blocks[order[i]] + by_boundary[block if block != self.NO_BLOCK_INDEX else i + 1][species] = gap - if by_boundary: - busiest = max(by_boundary.values(), key=len) - boundary_n = len(busiest) - spread = max(busiest.values()) - min(busiest.values()) - else: - boundary_n = 0 - spread = 0 + if by_boundary: + busiest = max(by_boundary.values(), key=len) + boundary_n = len(busiest) + spread = max(busiest.values()) - min(busiest.values()) + else: + boundary_n = 0 + spread = 0 - return [ - len(per_species_chunks), - counts["single"], counts["contiguous"], counts["split"], - counts["multi_scaffold"], counts["multi_strand"], counts["no_bases"], - max(per_species_total) if per_species_total else 0, - int(statistics.median(all_gaps)) if all_gaps else 0, - len(by_boundary), - boundary_n, - spread, - ] + return [ + len(self.species), + counts["single"], counts["contiguous"], counts["split"], + counts["multi_scaffold"], counts["multi_strand"], counts["no_bases"], + max(per_species_total) if per_species_total else 0, + int(statistics.median(all_gaps)) if all_gaps else 0, + len(by_boundary), + boundary_n, + spread, + ] + + +def summarizeElementLoci(per_species_chunks): + """ + ElementLociAccumulator over {species: [chunk dict, ...]}. Every key counts + toward n.species, even with an empty list. + """ + acc = ElementLociAccumulator() + for species, chunks in per_species_chunks.items(): + acc.addSpecies(species) + for c in chunks: + if c['size'] > 0: + acc.add(species, c['scaffold'], c['strand'], c['start'], c['size'], c.get('block_index')) + else: + acc.add(species, None, None, None, c['size'], None) + return acc.summary() def fastaScaffoldField(chunks): @@ -1399,8 +1510,13 @@ def fetchByRegion( block_lengths = [] # Holds lengths of each block written block_stats = [] # Holds [new_ref_start, new_ref_end, new_ref_size] for each block - # species -> [chunk, ...] for this element, feeding summarizeElementLoci. - element_chunks = defaultdict(list) + # Per-element locus summary. MAF output streams: each chunk is folded in as + # it arrives, so memory does not grow with the region (see + # ElementLociAccumulator). FASTA output must keep every chunk dict anyway + # (headers, --loci-table), so it only collects references to them and + # summarizes once at the end, as before -- streaming would save nothing. + element_loci = None if as_fasta else ElementLociAccumulator() + element_chunks = defaultdict(list) if as_fasta else None block_ordinal = 0 total_ref_bases = 0 # Total ref bases across all blocks @@ -1480,11 +1596,12 @@ def fetchByRegion( # while the FASTA path counts only the kept one. block_ordinal += 1 if not as_fasta: - # Mutated in place rather than copied: trimMafBlock built these - # dicts fresh for this block and nothing else reads them. + # trimMafBlock built these dicts for this block only; nothing keeps + # them, so they are freed with the block. + add_chunk = element_loci.add for chunk in block_chunks: - chunk['block_index'] = block_ordinal - element_chunks[chunk.pop('species')].append(chunk) + add_chunk(chunk['species'], chunk['scaffold'], chunk['strand'], + chunk['start'], chunk['size'], block_ordinal) if not single_output and trimmed is not None: block_stats.append(trimmed_info) @@ -1637,7 +1754,8 @@ def fetchByRegion( profile_state["regions"] += 1 addProfile(profile_state, "time_index_scan", time.perf_counter() - scan_timer_start) addProfile(profile_state, "time_region_total", time.perf_counter() - region_timer_start) - element_cols = "\t".join(str(v) for v in summarizeElementLoci(element_chunks)) + element_values = summarizeElementLoci(element_chunks) if as_fasta else element_loci.summary() + element_cols = "\t".join(str(v) for v in element_values) summary = (f"{region_str}\t{num_blocks}\t{blocks_length}\t{len(species_order) if as_fasta else 'NA'}\t{element_cols}") return summary diff --git a/tests/test_loci_accumulator.py b/tests/test_loci_accumulator.py new file mode 100644 index 0000000..1a6232c --- /dev/null +++ b/tests/test_loci_accumulator.py @@ -0,0 +1,260 @@ +""" +ElementLociAccumulator must compute exactly what the v0.8.1 per-element summary +computed from every chunk held in memory. The oracle below is that code, +copied verbatim from mafutils v0.8.1 (renamed oracle_*), so this test pins the +new streaming version to the old results value for value. +""" +import random +import statistics +from collections import defaultdict + +import pytest + +from mafutils.fetch import ELEMENT_LOCI_HEADERS, ElementLociAccumulator, summarizeElementLoci + + +# --------------------------------------------------------------------------- +# Oracle: verbatim from mafutils v0.8.1 fetch.py +# --------------------------------------------------------------------------- + +def oracle_classifyChunks(chunks): + """ + Classifies what a species' contribution to one region actually is, from its + per-block chunks. Returns (class, max_gap). + + A MAF's blocks are contiguous in the REFERENCE only. A non-reference + species' chunks may be separated, inverted, or on different source + sequences, and stitching concatenates them without padding -- so a stitched + row is not necessarily one locus in that genome. Measured on real data + (2,000 regions, 15 species): 77.8% single, 10.7% contiguous, 10.6% split, + 0.05% multi-strand, 0.8% multi-scaffold. + + max_gap is the largest-magnitude gap between consecutive chunks (sorted by + start), reported SIGNED -- negative means chunks overlap in the source, + which indicates something different again (paralogy, tandem repeat). It is + 0 for single/contiguous and None for multi_scaffold, where a gap between + positions on different sequences is meaningless. + + A split gap is NOT an indel the alignment absorbed. The aligner threads + query insertions inline as reference-row gap columns, but only up to ~37bp + on real data (156,928 such runs observed, max 37, none >=40); past that it + breaks the block instead. So a gap between chunks is query sequence the + aligner *declined to align* -- real bases absent from the emitted FASTA + while the header's bounding-box span still covers them. Consequence, + verified with zero exceptions across 27,599 records: single/contiguous have + header span == emitted bases, and split/multi_* never do. + + No gap threshold is applied on purpose: "how large a gap disqualifies a + locus" is a judgment for the analysis (real median split 45bp, max 335Mb), + so the number is reported and the caller decides. + + Chunks contributing zero aligned bases are excluded from the gap/class + arithmetic: a species can have an s-line in a block whose aligned sequence + falls entirely outside the requested region once trimmed, leaving an + all-gap row. If NO chunk contributes bases the class is "no_bases" rather + than "single" -- 5.46% of real rows are zero-size, and 1,648 of them used + to report as a clean "single" with a zero-width header span. + """ + if not chunks: + return "none", None + + contributing = [c for c in chunks if c['size'] > 0] + if not contributing: + return "no_bases", 0 + if len(contributing) == 1: + return "single", 0 + if len({c['scaffold'] for c in contributing}) > 1: + return "multi_scaffold", None + if len({c['strand'] for c in contributing}) > 1: + return "multi_strand", None + + ordered = sorted(contributing, key=lambda c: c['start']) + gaps = [ + ordered[i + 1]['start'] - (ordered[i]['start'] + ordered[i]['size']) + for i in range(len(ordered) - 1) + ] + if all(g == 0 for g in gaps): + return "contiguous", 0 + return "split", max(gaps, key=abs) + + +def oracle_chunkGaps(chunks): + """ + Signed gaps between a species' consecutive contributing chunks, in + ascending-start order. Shared by classifyChunks' caller and the per-element + summary so both see the same definition. + """ + contributing = sorted((c for c in chunks if c['size'] > 0), key=lambda c: c['start']) + if len(contributing) < 2: + return [] + if len({c['scaffold'] for c in contributing}) > 1 or len({c['strand'] for c in contributing}) > 1: + return [] + return [ + contributing[i + 1]['start'] - (contributing[i]['start'] + contributing[i]['size']) + for i in range(len(contributing) - 1) + ] + + +def oracle_summarizeElementLoci(per_species_chunks): + """ + Per-ELEMENT summary of how each species' contribution looks, as the values + for ELEMENT_LOCI_HEADERS. + + Classification is irreducibly per-species -- one real element was + contiguous for 5 species, split for 6 and empty for 3 -- so this reports + counts per class plus the extent and the concentration of the splits, + leaving every threshold to the caller. + + Extent needs two numbers, not one. A real element had a median gap of 63bp + across 13 split species but one species off by 335,617,958bp: the max alone + condemns an element where 12 species are nearly fine, the median alone + waves through a catastrophic misalignment. + - split.bases.max : max over species of total unaligned bases, i.e. + exactly (header span - aligned bases) for the worst + species -- how much sequence its row is missing. + - split.gap.median : median of every individual gap in the element. + + Concentration answers "is one shared event responsible?". Splits occur at + block boundaries, which are single reference positions, so species that + jump at the SAME boundary are jumping at the same place. Comparing + split.max.boundary.n.species to n.split distinguishes one ancestral indel + (equal, with a small gap spread -- a real 33/34/33/33/33 case gives spread + 1) from independent lineage-specific events (much lower). + + A species absent from a block inside its own span (1.6% of real multi-chunk + cases) has a gap spanning more than one boundary; it is attributed to the + first boundary it crosses, and the per-chunk loci table carries + block_index/gap_to_next so that case stays inspectable. + """ + counts = {"single": 0, "contiguous": 0, "split": 0, + "multi_scaffold": 0, "multi_strand": 0, "no_bases": 0} + all_gaps = [] + per_species_total = [] + # boundary -> {species: gap}, keyed by the block ordinal the gap starts at + by_boundary = defaultdict(dict) + + for species, chunks in per_species_chunks.items(): + locus_class, _max_gap = oracle_classifyChunks(chunks) + if locus_class in counts: + counts[locus_class] += 1 + + gaps = oracle_chunkGaps(chunks) + nonzero = [g for g in gaps if g != 0] + if nonzero: + all_gaps.extend(nonzero) + per_species_total.append(sum(nonzero)) + contributing = sorted((c for c in chunks if c['size'] > 0), key=lambda c: c['start']) + for i, gap in enumerate(gaps): + if gap != 0: + by_boundary[contributing[i].get('block_index', i + 1)][species] = gap + + if by_boundary: + busiest = max(by_boundary.values(), key=len) + boundary_n = len(busiest) + spread = max(busiest.values()) - min(busiest.values()) + else: + boundary_n = 0 + spread = 0 + + return [ + len(per_species_chunks), + counts["single"], counts["contiguous"], counts["split"], + counts["multi_scaffold"], counts["multi_strand"], counts["no_bases"], + max(per_species_total) if per_species_total else 0, + int(statistics.median(all_gaps)) if all_gaps else 0, + len(by_boundary), + boundary_n, + spread, + ] + + +# --------------------------------------------------------------------------- +# Random elements that hit every branch +# --------------------------------------------------------------------------- + + +def random_element(rng): + """ + {species: [chunk, ...]} with start ties, out-of-order and overlapping + chunks, zero sizes, several scaffolds and strands (in either order), + repeated and missing block_index, and species with no chunks at all. + """ + element = {} + for s in range(rng.randint(1, 8)): + chunks = [] + kind = rng.random() + scaffolds = ["s1"] if kind < 0.6 else ["s1", "s2"] + strands = ["+"] if kind < 0.4 or kind > 0.8 else ["+", "-"] + pos = rng.randint(0, 50) + for _ in range(rng.randint(0, 10)): + if rng.random() < 0.5: # mostly collinear... + start = pos + rng.choice([0, 0, 0, rng.randint(-5, 40)]) + else: # ...but not always + start = rng.randint(0, 200) + size = 0 if rng.random() < 0.15 else rng.randint(1, 30) + pos = start + size + c = {"scaffold": rng.choice(scaffolds), "start": start, "size": size, + "end": start + size, "strand": rng.choice(strands), "src_size": 1000} + if rng.random() > 0.05: + c["block_index"] = rng.randint(1, 6) + chunks.append(c) + element[f"sp{s}"] = chunks + return element + + +@pytest.mark.parametrize("seed", range(20)) +def test_matches_v081_on_random_elements(seed): + rng = random.Random(seed) + for _ in range(500): + element = random_element(rng) + expected = oracle_summarizeElementLoci(element) + assert summarizeElementLoci(element) == expected, element + + +def test_matches_v081_when_fed_chunk_by_chunk_in_block_order(): + """fetch adds chunks block by block, interleaving species -- not per species.""" + rng = random.Random(1234) + for _ in range(2000): + element = random_element(rng) + events = [(c.get("block_index", 0), sp, i, c) + for sp, chunks in element.items() for i, c in enumerate(chunks)] + events.sort(key=lambda e: (e[0], e[2])) + # the order fetch would see them in, and the per-species dicts in that order + acc = ElementLociAccumulator() + reordered = defaultdict(list) + for _block, sp, _i, c in events: + reordered[sp].append(c) + if c["size"] > 0: + acc.add(sp, c["scaffold"], c["strand"], c["start"], c["size"], c.get("block_index")) + else: + acc.add(sp, None, None, None, c["size"], None) + assert acc.summary() == oracle_summarizeElementLoci(dict(reordered)), dict(reordered) + + +def test_multi_species_drop_their_chunk_arrays(): + """The memory fix: no per-chunk storage once gaps can no longer matter.""" + acc = ElementLociAccumulator() + for i in range(1000): + acc.add("collinear", "s1", "+", i * 10, 10, i) + acc.add("rearranged", "s1" if i < 500 else "s2", "+", i * 10, 10, i) + acc.add("inverted", "s1", "+" if i < 500 else "-", i * 10, 10, i) + assert len(acc.species["collinear"][6]) == 1000 + assert acc.species["rearranged"][6] is None + assert acc.species["inverted"][6] is None + vals = dict(zip(ELEMENT_LOCI_HEADERS, acc.summary())) + assert (vals["n.contiguous"], vals["n.multi.scaffold"], vals["n.multi.strand"]) == (1, 1, 1) + + +def test_scaffold_change_after_strand_change_is_multi_scaffold(): + """multi_scaffold outranks multi_strand even after the arrays are gone.""" + acc = ElementLociAccumulator() + acc.add("sp", "s1", "+", 0, 10, 1) + acc.add("sp", "s1", "-", 20, 10, 2) # multi_strand: arrays dropped + acc.add("sp", "s2", "-", 40, 10, 3) # ...then a second scaffold + element = {"sp": [ + {"scaffold": "s1", "strand": "+", "start": 0, "size": 10, "block_index": 1}, + {"scaffold": "s1", "strand": "-", "start": 20, "size": 10, "block_index": 2}, + {"scaffold": "s2", "strand": "-", "start": 40, "size": 10, "block_index": 3}, + ]} + assert acc.summary() == oracle_summarizeElementLoci(element) + assert dict(zip(ELEMENT_LOCI_HEADERS, acc.summary()))["n.multi.scaffold"] == 1 From bc9c65dd2c914572afdf80a5bd3115660ebbd4a0 Mon Sep 17 00:00:00 2001 From: gwct Date: Sun, 8 Nov 2026 21:08:45 -0500 Subject: [PATCH 2/3] Count only s lines as sequences in the index (index format 3) Since the first release, the index's per-block sequence count (block.idx column 6, and the seq_lines/max_seqs header totals) included i, e and q lines. Cactus MAFs have none of these lines, but UCSC multiz MAFs do. On a 45-way multiz MAF, v0.8.1 reported max_seqs=89 and 633M sequence lines; the true values are 45 and 270M. index: - Counts only s lines, reference included. - Stops with one error, and writes no index, on an unrecognized line type, an s line before the first a line, or a block with no s line. - 5-10% faster on every input measured (84 GB multiz MAF: 323 -> 290 s). Memory is unchanged. - Gives the same output otherwise: on a 42.7 GB cactus MAF (uncompressed, gz, bgzip) the index bodies are byte-identical to v0.8.1. On the multiz MAF only column 6 changes, and it matches an independent s-line count for all 8.4M blocks. Index format 3: - New indexes say format=3, because column 6 changed meaning. - stats (unless --no-block-table) and info refuse an index of format 2 or lower, or one with no header, and say how to rebuild it. validate reports the format without changing its verdict. fetch and gc never read the count, so they still accept old indexes. Tests: - New tests/test_index_seq_lines.py, on tests/multiz45way-excerpt.maf (14 real blocks with 72 s, 58 i and 12 e lines). - The 13 committed test indexes were regenerated, with bodies confirmed identical, and their headers changed from format=2 to format=3. test_validate and test_info now expect format 3. Co-Authored-By: Claude Opus 5.5 --- DEVELOPMENT.md | 78 +++++- README.md | 9 +- mafutils/index.py | 30 ++- mafutils/info.py | 20 +- mafutils/lib/common.py | 157 ++++++++--- mafutils/stats.py | 8 + mafutils/validate.py | 8 + tests/example-tabs.maf.block.idx | 2 +- tests/example-tabs.maf.scaffold.idx | 2 +- tests/example.maf.bgz.block.idx | 2 +- tests/example.maf.bgz.scaffold.idx | 2 +- tests/example.maf.block.idx | 2 +- tests/example.maf.gz.block.idx | 2 +- tests/example.maf.gz.scaffold.idx | 2 +- tests/example.maf.scaffold.idx | 2 +- tests/example.maf.scaffold.regenerated.idx | 2 +- tests/loci-classes.maf.block.idx | 2 +- tests/loci-classes.maf.scaffold.idx | 2 +- tests/multiz45way-excerpt.maf | 171 ++++++++++++ tests/real-excerpt.maf.block.idx | 2 +- tests/real-excerpt.maf.scaffold.idx | 2 +- tests/test_index_seq_lines.py | 294 +++++++++++++++++++++ tests/test_info.py | 9 +- tests/test_validate.py | 2 +- 24 files changed, 731 insertions(+), 81 deletions(-) create mode 100644 tests/multiz45way-excerpt.maf create mode 100644 tests/test_index_seq_lines.py diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index 5cf045e..74fd1b6 100644 --- a/DEVELOPMENT.md +++ b/DEVELOPMENT.md @@ -35,9 +35,10 @@ for user-facing installation/usage docs. 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"`. + MISMATCH -- the same dict goes into both. The keys did not change the format + (it stayed at 2): they are purely additive, and `readIndexHeader` already + parses arbitrary keys. Format 3 came later, for a change in what the + sequence count *means* (see "count only `s` lines" below). `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 @@ -51,6 +52,58 @@ 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. +- **The index's per-block sequence count is `s` lines only (index format + 3).** + - **The bug:** since the first commit, `iterMafBlocks` kept every line that + was not `a`, `#` or blank, and `num_seqs` was `len(block) - 1`. That + counted `i`, `e` and `q` lines as sequences in block.idx column 6 and in + the `seq_lines`/`max_seqs` totals. Cactus MAFs have none of these lines, + so this went unnoticed. On a 45-way UCSC-style multiz MAF (84 GB, 8.4M + blocks), v0.8.1 reported `seq_lines=633,354,590 max_seqs=89`; the true + values are 270,212,333 and 45. 7.9M blocks claimed more than 45 sequences. + phyloacc-workflows chunks on this column, so its chunks were far emptier + than intended. + - **The fix:** `iterMafBlocks` tests for `s` first and only counts the + non-reference ones, yielding `(ref_line, n_seq_lines)`. `i`/`e`/`q` lines + are skipped by one tuple `startswith`. Any other line type, an `s` line + before the first `a`, or a block with no `s` line raises `MafFormatError`. + `index` reports it as one error and writes no index (the final files are + only written after the scan). + - **Why the format changed:** the column's *meaning* changed, which no + reader can detect from the values. `stats` (when it writes `block.tsv`, + whose `n_seq_lines_index` is this count) and `info` refuse a format-2 or + headerless index. `validate` reports it without changing its verdict. + `fetch` and `gc` never read the count, so they accept old indexes. This + avoids forcing a rebuild of, e.g., a 7.4 TB MAF's index just to fetch + from it. + - **Faster:** the old loop paid a `startswith('#')`, a `strip()` that copied + the whole line, a `startswith('a')` and a list append for every sequence + line. Now an `s` line costs one `startswith` and an increment. + `mafutils index`, v0.8.1 vs new (single process, `/usr/bin/time`, pairs + with the run order alternated): + + | Input | v0.8.1 | New | Change | + |---|---|---|---| + | multiz 84 GB, uncompressed, 3 pairs | 323.4s | 289.8s | -10% | + | hamster 42.7 GB, uncompressed, 2 pairs | 123.9s | 110.9s | -10% | + | hamster gz, 1 pair | 184.3s | 169.8s | -8% | + | hamster bgzip, 1 pair | 226.1s | 215.2s | -5% | + | 1.5 GB slice gz, 3 pairs | 7.1s | 6.5s | -8% | + | 1.5 GB slice bgzip, 3 pairs | 8.4s | 7.8s | -7% | + + The ranges of the two versions never overlapped. Peak RSS stayed at + 31-46 MB. + - **Same output otherwise:** + - hamster (none/gz/bgzip) and the slices: block and scaffold index bodies + are byte-identical to v0.8.1, and the headers differ only in `format=3`; + - multiz: every column but 6 and the scaffold index are identical, and + column 6 equals an independent `awk` count of `s` lines for all + 8,436,088 blocks. + - **Tests:** `tests/test_index_seq_lines.py`, on + `tests/multiz45way-excerpt.maf` (14 real blocks with 72 `s`, 58 `i` and + 12 `e` lines). The committed test indexes were regenerated, their bodies + confirmed identical (those MAFs have only `s` lines), and their headers + changed to `format=3`. - **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` @@ -230,15 +283,18 @@ for user-facing installation/usage docs. before assuming it's mafutils-specific. Keep this file otherwise unmodified so future diffs against upstream stay clean. - Index files (`.block.idx`/`.scaffold.idx`) carry a `# mafutils-index - format=2 maf= compression= size= mtime= + format=3 maf= compression= size= mtime= hash=:` header line (helpers in `mafutils/lib/common.py`: `writeIndexHeader`/`readIndexHeader`/`compareIndexHeader`/ - `validateIndexHeader`). Existing parsers already skip `#` lines, so this - is backward compatible. Format-1 indexes (pre-dating size/mtime/hash) - simply lack those keys in the parsed header dict — every consumer treats - their absence as "can't check this dimension," not an error, same as a - fully missing header. Tests build a headerless index on the fly rather than - relying on a committed one (see above). + `validateIndexHeader`; format history at `INDEX_FORMAT`). Existing parsers + already skip `#` lines, so this is backward compatible. Format-1 indexes + (pre-dating size/mtime/hash) simply lack those keys in the parsed header + dict — every consumer treats their absence as "can't check this + dimension," not an error, same as a fully missing header. The exception is + the per-block sequence count, which formats 1-2 got wrong on MAFs with + `i`/`e`/`q` lines: `stats` (block table) and `info` refuse those indexes. + Tests build a headerless index on the fly rather than relying on a + committed one (see above). - **Computing the hash costs nothing extra at `mafutils index` time, for all three compression types** — this was the whole point of always storing it. `common._HashingRawIO` (an `io.RawIOBase` subclass) sits @@ -343,7 +399,7 @@ for user-facing installation/usage docs. C-implemented `BufferedReader.readline()` while adding Python-level buffer slicing), and a block-level scanner that finds `b"\na"` boundaries to avoid per-line iteration entirely (2.71s — still has to walk every line to count - non-blank lines for `num_seqs`, plus buffer-concatenation cost). Remaining + the `s` lines for `num_seqs`, plus buffer-concatenation cost). Remaining headroom is real but modest: end-to-end sits at ~129 MB/s against a ~640 MB/s ceiling (raw read + MD5 with no parsing at all), and closing it would require assuming things about intra-block blank/comment lines that risk silently diff --git a/README.md b/README.md index 5db5f97..29c392a 100644 --- a/README.md +++ b/README.md @@ -143,6 +143,13 @@ 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). +**Older indexes.** An index whose header says `format=2` or lower, or that +has no header, counted `i`/`e`/`q` lines as sequences. That inflates its +per-block sequence count on MAFs that have such lines (UCSC multiz MAFs do; +cactus MAFs don't). `fetch` and `gc` don't use the count and still accept +these indexes. `stats` (unless `--no-block-table`) and `info` show it, so they +ask you to rebuild with `mafutils index`, and `validate` flags it. + ## Commands ### `mafutils index` @@ -222,7 +229,7 @@ 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 -reference. `aln columns` is the block width, gap columns included, which all +reference (`i`/`e`/`q` lines are not sequences). `aln columns` is the block width, gap columns included, which all species share. ### `mafutils fetch` diff --git a/mafutils/index.py b/mafutils/index.py index d92962d..271fe39 100644 --- a/mafutils/index.py +++ b/mafutils/index.py @@ -9,6 +9,7 @@ import logging import os import shutil +import sys import tempfile from typing import Annotated, Optional @@ -19,23 +20,25 @@ from mafutils.lib import options as OPTIONS -def process_maf_block(block): +def process_maf_block(ref_line, num_seqs): """ - Returns (row, ref_len, aln_cols, num_seqs) for one block: the index row as - strings, plus the three integers the header aggregates need. + Returns (row, ref_len, aln_cols) for one block, from its reference 's' + line and its count of 's' lines (both from iterMafBlocks): the index row + as strings, plus the two 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_seq = ref_line.split() ref_scaff = ref_seq[1].split(".", 1)[1] - line_len = str(len(block[1])) - num_seqs = len(block) - 1 # every s line, including the reference + line_len = str(len(ref_line)) seq_len = len(ref_seq[6]) # block width, gap columns included + # num_seqs: 's' lines only, the reference included. i/e/q lines are not + # sequences (index format 3; see COMMON.INDEX_FORMAT). 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 + return row, int(ref_seq[3]), seq_len def run_index(maf_file, block_index_path, scaffold_index_path, LOG, progress_interval=COMMON.PROGRESS_INTERVAL_SECONDS): @@ -82,8 +85,8 @@ def run_index(maf_file, block_index_path, scaffold_index_path, LOG, progress_int max_seqs = 0 distinct_scaffolds = set() - for block, block_start, block_end in COMMON.iterMafBlocks(maf_stream, binary=use_binary): - block_info, ref_len, aln_cols, num_seqs = process_maf_block(block) + for (ref_line, num_seqs), block_start, block_end in COMMON.iterMafBlocks(maf_stream, binary=use_binary): + block_info, ref_len, aln_cols = process_maf_block(ref_line, num_seqs) ref_scaffold = block_info[0] n_blocks += 1 @@ -163,4 +166,11 @@ def index_command( log_verbosity="SCREEN", logger_name="maf_index_logger", ) - run_index(maf_file, block_index, scaffold_index, logging.getLogger(logger_name), progress_interval=progress_interval) + LOG = logging.getLogger(logger_name) + try: + run_index(maf_file, block_index, scaffold_index, LOG, progress_interval=progress_interval) + except COMMON.MafFormatError as e: + # Raised mid-scan, before run_index writes the final index files, so + # nothing is left behind for a MAF that failed. + LOG.error(f"Cannot index {maf_file}: {e}. No index was written.") + sys.exit(1) diff --git a/mafutils/info.py b/mafutils/info.py index 3df5694..9377c61 100644 --- a/mafutils/info.py +++ b/mafutils/info.py @@ -21,6 +21,10 @@ recommends rebuilding, so a long wait is never a surprise. 3. No index: report only what needs no index and point at `mafutils index`. +An index older than format 3 (see COMMON.INDEX_FORMAT) may have counted i/e/q +lines as sequences, so its `seq lines` / `max seqs/block` are shown as +unavailable and info exits 1 after printing everything else. + 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 @@ -144,6 +148,7 @@ def run_info(args, LOG): header = COMMON.readIndexHeader(index_file) if have_index else None aggregates = readAggregates(header) + seq_counts_exact = not have_index or COMMON.seqCountsExact(header) if have_index and aggregates is None: index_size = os.path.getsize(index_file) @@ -201,8 +206,14 @@ def run_info(args, LOG): 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']:,}") + if seq_counts_exact: + 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: + old_index = f"index format {COMMON.indexFormat(header)}" if header else "a headerless index" + unavailable = f"unavailable ({old_index} may count i/e/q lines)" + out.append(f"{'seq lines':<15}: {unavailable}") + out.append(f"{'max seqs/block':<15}: {unavailable}") else: out.append(f"{'blocks':<15}: unavailable without an index -- run `mafutils index` first") @@ -227,6 +238,11 @@ def run_info(args, LOG): print("\n".join(out)) + if not seq_counts_exact: + sys.stdout.flush() # the report first, then the error under it + LOG.error(COMMON.oldIndexSeqCountMessage(index_file, args.maf_file, header)) + sys.exit(1) + ############################################################################# diff --git a/mafutils/lib/common.py b/mafutils/lib/common.py index 8dec4c2..48b9d38 100644 --- a/mafutils/lib/common.py +++ b/mafutils/lib/common.py @@ -307,6 +307,20 @@ def deriveScaffoldIndexPath(maf_file): # them simply lacks the keys and consumers fall back to scanning. INDEX_AGGREGATE_KEYS = ("blocks", "scaffolds", "ref_bases", "aln_cols", "seq_lines", "max_seqs") +# Format history: +# 1 no size/mtime/hash in the header +# 2 size/mtime/hash (and, from v0.7.0, the aggregate keys -- additive, so +# they did not need a new format) +# 3 the per-block sequence count (block.idx column 6) and the seq_lines / +# max_seqs totals count only 's' lines. Formats 1-2 also counted i/e/q +# lines, so on a MAF that has them (UCSC multiz MAFs do; cactus MAFs do +# not) those numbers are inflated -- up to ~2x on a 45-way multiz MAF. +# Block offsets and every other column are unchanged. +INDEX_FORMAT = 3 + +# The first format whose sequence counts are 's' lines only. +SEQ_COUNT_INDEX_FORMAT = 3 + def writeIndexHeader(stream, maf_file, compression, size, mtime, content_hash, aggregates=None): """ @@ -319,10 +333,10 @@ def writeIndexHeader(stream, maf_file, compression, size, mtime, content_hash, a 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 key=value tokens are additive (readIndexHeader parses arbitrary keys, + and every consumer treats a missing key as "can't check this dimension"), + so they do not change the format. A change in what an existing column + means does -- see INDEX_FORMAT. """ extra = "" if aggregates: @@ -330,11 +344,42 @@ def writeIndexHeader(stream, maf_file, compression, size, mtime, content_hash, a 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"{INDEX_HEADER_PREFIX} format={INDEX_FORMAT} maf={os.path.basename(maf_file)} " f"compression={compression} size={size} mtime={mtime} hash={content_hash}{extra}\n" ) +def indexFormat(header): + """The index's format number from a readIndexHeader dict; 0 for a missing header.""" + if header is None: + return 0 + value = header.get("format") + if value is None or not value.isdigit(): + raise ValueError(f"index header has no valid format number: format={value}") + return int(value) + + +def seqCountsExact(header): + """True when the index's sequence counts are 's' lines only (format 3+).""" + return indexFormat(header) >= SEQ_COUNT_INDEX_FORMAT + + +def oldIndexSeqCountMessage(index_file, maf_file, header): + """ + The one wording stats/info/validate use for an index whose sequence counts + predate format 3. Says "may": the counts are only wrong if the MAF has + i/e/q lines, which the index alone cannot tell. + """ + built_by = ( + "an older mafutils with no index header" if header is None + else f"an older mafutils (index format {indexFormat(header)})" + ) + return ( + f"{index_file} was built by {built_by}, whose per-block sequence counts " + f"may also include i/e/q lines. Rebuild it with `mafutils index {maf_file}`." + ) + + def readIndexHeader(index_file): """ Peeks at the first line of an index file. Returns a dict of the header's @@ -498,23 +543,42 @@ def validateIndexHeader(header, maf_file, detected_compression, LOG, strict=Fals ############################################################################# +class MafFormatError(ValueError): + """A MAF line or block that mafutils cannot index correctly.""" + + +def _linePreview(line): + if isinstance(line, bytes): + line = line.decode("utf-8", errors="replace") + line = line.rstrip("\r\n") + return repr(line[:60] + ("..." if len(line) > 60 else "")) + + def iterMafBlocks(stream, binary=False): """ - Splits an open MAF stream into raw per-block lines, split on 'a'-prefixed - header lines (skipping blank/comment lines). Yields - (block_lines, block_start_offset, block_end_offset) tuples, where the - offsets bracket each block and are suitable for later random-access - seeking back into the same file. - - block_lines[0] (the 'a' header) and block_lines[1] (the reference/first - sequence line) are always returned as stripped `str` -- they are the - only two lines whose *content* mafutils index reads. The remaining - sequence lines are returned exactly as readline() produced them, and in - binary mode are still `bytes`: only their *count* (len(block_lines)) is - ever used downstream, so decoding/stripping them would be pure overhead. - That matters enormously at scale -- a block can carry hundreds of - sequence lines (one per species) and a whole-genome MAF hundreds of - millions of blocks. + Splits an open MAF stream into blocks at 'a' lines. Yields + ((ref_line, n_seq_lines), block_start_offset, block_end_offset), where + the offsets bracket each block (any comment/blank lines after it + included) for later random-access seeking back into the same file. + + ref_line is the block's first 's' line (the reference), stripped `str`: + the only line whose *content* mafutils index reads. n_seq_lines counts + the block's 's' lines, the reference included. The other 's' lines are + only counted -- never stored, decoded or stripped -- because a block can + carry hundreds of them and a whole-genome MAF has hundreds of millions + of blocks. + + Line types, tested in order of frequency so the common 's' line costs a + single startswith(): + s counted (the first one in a block is the reference) + a starts a new block + i, e, q skipped: per-species context, empty-region and quality lines + describe a species' sequence but are not one. Counting them + (as mafutils did before index format 3) inflated the per-block + count -- on a UCSC 45-way multiz MAF, to 89 "sequences". + #, blank skipped + Anything else, an 's' line before the first 'a' line, or a block with no + 's' line raises MafFormatError rather than being guessed at. Two position-tracking modes, because the two cases are fundamentally different: @@ -540,14 +604,17 @@ def iterMafBlocks(stream, binary=False): why the binary counting path above must not be used for bgzip. """ if binary: - empty, hash_prefix, a_prefix = b"", b"#", b"a" + empty, hash_prefix, a_prefix, s_prefix = b"", b"#", b"a", b"s" + skipped_prefixes = (b"i", b"e", b"q") else: - empty, hash_prefix, a_prefix = "", "#", "a" + empty, hash_prefix, a_prefix, s_prefix = "", "#", "a", "s" + skipped_prefixes = ("i", "e", "q") - block_lines = [] + ref_line = None + n_seq_lines = 0 + n_blocks = 0 block_start = None pos = 0 if binary else stream.tell() - next_line_is_reference = False while True: line = stream.readline() @@ -555,26 +622,38 @@ def iterMafBlocks(stream, binary=False): break next_pos = pos + len(line) if binary else stream.tell() - if line.startswith(hash_prefix) or not line.strip(): - pos = next_pos - continue - - if line.startswith(a_prefix): - if block_lines: - yield block_lines, block_start, pos + if line.startswith(s_prefix): + if n_seq_lines: + n_seq_lines += 1 + else: + if block_start is None: + raise MafFormatError(f"'s' line before the first 'a' line: {_linePreview(line)}") + ref_line = line.decode().strip() if binary else line.strip() + n_seq_lines = 1 + elif line.startswith(a_prefix): + if block_start is not None: + if not n_seq_lines: + raise MafFormatError(f"block {n_blocks:,} has no 's' line") + yield (ref_line, n_seq_lines), block_start, pos + n_blocks += 1 block_start = pos - block_lines = [line.decode().strip() if binary else line.strip()] - next_line_is_reference = True - elif next_line_is_reference: - block_lines.append(line.decode().strip() if binary else line.strip()) - next_line_is_reference = False + ref_line = None + n_seq_lines = 0 + elif line.startswith(skipped_prefixes) or line.startswith(hash_prefix) or not line.strip(): + pass else: - block_lines.append(line) + where = f"in block {n_blocks:,}" if n_blocks else "before the first block" + raise MafFormatError( + f"unrecognized line {where}: {_linePreview(line)} " + f"(expected a, s, i, e, q, # or blank lines)" + ) pos = next_pos - if block_lines: - yield block_lines, block_start, pos if binary else stream.tell() + if block_start is not None: + if not n_seq_lines: + raise MafFormatError(f"block {n_blocks:,} has no 's' line") + yield (ref_line, n_seq_lines), block_start, pos if binary else stream.tell() ############################################################################# diff --git a/mafutils/stats.py b/mafutils/stats.py index 77bdec6..51c9821 100644 --- a/mafutils/stats.py +++ b/mafutils/stats.py @@ -1492,6 +1492,14 @@ def run_stats(args, cmdline="mafutils stats"): LOG.info(f"Reading index file: {index_file}") index_header = COMMON.readIndexHeader(index_file) COMMON.validateIndexHeader(index_header, args.maf_file, maf_compression, LOG, strict=args.verify_hash) + # block.tsv's n_seq_lines_index column is the index's per-block count, the + # only stats output that uses it (the dashboard and overall.tsv do not). + if not args.no_block_table and not COMMON.seqCountsExact(index_header): + LOG.error( + f"{COMMON.oldIndexSeqCountMessage(index_file, args.maf_file, index_header)} " + f"block.tsv reports these counts (n_seq_lines_index); --no-block-table runs without them." + ) + sys.exit(1) header_blocks = COMMON.headerBlockCount(index_header, LOG, progress=args.progress_interval > 0) LOG.info(COMMON.progressStartupLine(args.progress_interval)) diff --git a/mafutils/validate.py b/mafutils/validate.py index a231b54..f8e171d 100644 --- a/mafutils/validate.py +++ b/mafutils/validate.py @@ -93,6 +93,14 @@ def run_validate(args, cmdline="mafutils validate"): ) sys.exit(2) + # Reported, not part of the verdict: the verdict is whether the index + # matches this MAF, which an older format does not change. stats and info + # refuse it themselves where they would show its sequence counts. + if COMMON.seqCountsExact(header): + LOG.info(f"{'format':11}: {COMMON.indexFormat(header)}") + else: + LOG.warning(f"{'format':11}: {COMMON.oldIndexSeqCountMessage(index_file, args.maf_file, header)}") + result = COMMON.compareIndexHeader(header, args.maf_file, detected_compression, check_hash=True) conclusive_mismatch = False diff --git a/tests/example-tabs.maf.block.idx b/tests/example-tabs.maf.block.idx index b0c06c7..9a4b36f 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=778 mtime=1792526707.0 hash=md5:48230e6b5535b46c5e808c40e2dbe3f2 blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 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 diff --git a/tests/example-tabs.maf.scaffold.idx b/tests/example-tabs.maf.scaffold.idx index d5e7762..5de61c4 100644 --- a/tests/example-tabs.maf.scaffold.idx +++ b/tests/example-tabs.maf.scaffold.idx @@ -1,4 +1,4 @@ -# 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 +# mafutils-index format=3 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 diff --git a/tests/example.maf.bgz.block.idx b/tests/example.maf.bgz.block.idx index 17dbd08..96c3f8c 100644 --- a/tests/example.maf.bgz.block.idx +++ b/tests/example.maf.bgz.block.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example.maf.bgz compression=bgzip size=305 mtime=1792526737.0 hash=md5:54d36240911e394249124a0800fb75dd blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 maf=example.maf.bgz compression=bgzip size=305 mtime=1792526737.0 hash=md5:54d36240911e394249124a0800fb75dd 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 diff --git a/tests/example.maf.bgz.scaffold.idx b/tests/example.maf.bgz.scaffold.idx index 2374f60..84b54ca 100644 --- a/tests/example.maf.bgz.scaffold.idx +++ b/tests/example.maf.bgz.scaffold.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example.maf.bgz compression=bgzip size=305 mtime=1792526737.0 hash=md5:54d36240911e394249124a0800fb75dd blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 maf=example.maf.bgz compression=bgzip size=305 mtime=1792526737.0 hash=md5:54d36240911e394249124a0800fb75dd 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 diff --git a/tests/example.maf.block.idx b/tests/example.maf.block.idx index 68495d7..da1f2cc 100644 --- a/tests/example.maf.block.idx +++ b/tests/example.maf.block.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example.maf compression=none size=777 mtime=1792526707.0 hash=md5:5f78f4cd6b51323e127d2ab7f2e1e9d9 blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 maf=example.maf compression=none size=777 mtime=1792526707.0 hash=md5:5f78f4cd6b51323e127d2ab7f2e1e9d9 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 diff --git a/tests/example.maf.gz.block.idx b/tests/example.maf.gz.block.idx index fab89ef..b5cccf5 100644 --- a/tests/example.maf.gz.block.idx +++ b/tests/example.maf.gz.block.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example.maf.gz compression=gz size=281 mtime=1792526737.0 hash=md5:02ccf92d78b5dbcb6c80572b93bd4f9a blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 maf=example.maf.gz compression=gz size=281 mtime=1792526737.0 hash=md5:02ccf92d78b5dbcb6c80572b93bd4f9a 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 diff --git a/tests/example.maf.gz.scaffold.idx b/tests/example.maf.gz.scaffold.idx index 2ce8ea4..94a1f3b 100644 --- a/tests/example.maf.gz.scaffold.idx +++ b/tests/example.maf.gz.scaffold.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example.maf.gz compression=gz size=281 mtime=1792526737.0 hash=md5:02ccf92d78b5dbcb6c80572b93bd4f9a blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 maf=example.maf.gz compression=gz size=281 mtime=1792526737.0 hash=md5:02ccf92d78b5dbcb6c80572b93bd4f9a 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 diff --git a/tests/example.maf.scaffold.idx b/tests/example.maf.scaffold.idx index d18f436..3b93c7d 100644 --- a/tests/example.maf.scaffold.idx +++ b/tests/example.maf.scaffold.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example.maf compression=none size=777 mtime=1792526707.0 hash=md5:5f78f4cd6b51323e127d2ab7f2e1e9d9 blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 maf=example.maf compression=none size=777 mtime=1792526707.0 hash=md5:5f78f4cd6b51323e127d2ab7f2e1e9d9 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 diff --git a/tests/example.maf.scaffold.regenerated.idx b/tests/example.maf.scaffold.regenerated.idx index d18f436..3b93c7d 100644 --- a/tests/example.maf.scaffold.regenerated.idx +++ b/tests/example.maf.scaffold.regenerated.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=example.maf compression=none size=777 mtime=1792526707.0 hash=md5:5f78f4cd6b51323e127d2ab7f2e1e9d9 blocks=10 scaffolds=5 ref_bases=59 aln_cols=61 seq_lines=20 max_seqs=2 +# mafutils-index format=3 maf=example.maf compression=none size=777 mtime=1792526707.0 hash=md5:5f78f4cd6b51323e127d2ab7f2e1e9d9 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 diff --git a/tests/loci-classes.maf.block.idx b/tests/loci-classes.maf.block.idx index d1600fa..34e00a9 100644 --- a/tests/loci-classes.maf.block.idx +++ b/tests/loci-classes.maf.block.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=loci-classes.maf compression=none size=29422 mtime=1789746613.0 hash=md5:e2a9c32b01c9e0b4f798534f8815d6a2 +# mafutils-index format=3 maf=loci-classes.maf compression=none size=29422 mtime=1789746613.0 hash=md5:e2a9c32b01c9e0b4f798534f8815d6a2 CM000994.3 3541841 56 59 108 14 16 1666 CM000994.3 3541897 163 179 229 14 1666 5017 CM000994.3 3602158 290 311 361 15 5017 10577 diff --git a/tests/loci-classes.maf.scaffold.idx b/tests/loci-classes.maf.scaffold.idx index 3518fd5..f1bff58 100644 --- a/tests/loci-classes.maf.scaffold.idx +++ b/tests/loci-classes.maf.scaffold.idx @@ -1,2 +1,2 @@ -# mafutils-index format=2 maf=loci-classes.maf compression=none size=29422 mtime=1789746613.0 hash=md5:e2a9c32b01c9e0b4f798534f8815d6a2 +# mafutils-index format=3 maf=loci-classes.maf compression=none size=29422 mtime=1789746613.0 hash=md5:e2a9c32b01c9e0b4f798534f8815d6a2 CM000994.3 16 29422 diff --git a/tests/multiz45way-excerpt.maf b/tests/multiz45way-excerpt.maf new file mode 100644 index 0000000..8a1f3a4 --- /dev/null +++ b/tests/multiz45way-excerpt.maf @@ -0,0 +1,171 @@ +##maf version=1 scoring=roast +a score=19929.000000 +s galGal6.chr6 20883553 61 + 36374701 ----------------------ATTCTGTATAGCCGACAT----TACCTGCCTCGCATTCCTTGCCATtccc-----cccttcccacccccc +s HLtytAlb2.CABFVR010005803 1318225 49 + 2295574 ----------------------CTTTCGTACAGCTGGCACCGCTCGCTTACATT-CGTTCCTTGCCATCTCC-------------------- +i HLtytAlb2.CABFVR010005803 I 20 C 0 +s halLeu1.KL869536 4538483 45 - 5862070 ----------------------CTTTTGTATAGCCAGTAC----CGCTTACATT-CATTCCTTGCCCTCTCC-------------------- +i halLeu1.KL869536 I 20 C 0 +s aptFor1.KL225800 1006355 50 + 2472740 ----------------------CTTTTGTAGAGCCGGCAC----CCCTTACATT-CATTCTTTGCCATCTCCCCCTT--------------- +i aptFor1.KL225800 I 20 C 0 +s HLcalAnn5.CM012136 503863 44 + 6461084 -----------------------TTTGTTACAGCTGGCAC----CATTTATGTT-CACTCCTTGCCCTCTCC-------------------- +i HLcalAnn5.CM012136 C 0 C 0 +s cucCan1.KL448280 2754881 67 - 4414095 ttaaaaaaaagcaagcttctTGCTCTGGTGTAGCCAGCAC----TGCTTACATT-CATTCCTCGacatttac-------------------- +i cucCan1.KL448280 C 0 C 0 +e HLcalPug1.KQ482592 384241 439 + 1956262 I +e HLstrHab1.CM013776 2519986 0 - 7441144 C +e HLlorGal1.manual_scaffold_27 1988627 80 + 7065481 I +e HLfalTin1.ML997056 33510669 14 + 37096019 I +e falChe1.KB397934 1119933 13 + 4640569 I +e falPer1.KB390927 1430784 14 - 2055886 I + +a score=-11523.000000 +s galGal6.chr6 20883614 33 + 36374701 gccccccgcctttctcttccaaggtgacaaagg +s falChe1.KB397934 1119946 20 + 4640569 -----ACCCCTTGCTTTTGTGT--------AGC +i falChe1.KB397934 I 13 C 0 +s HLfalTin1.ML997056 33510683 18 + 37096019 -------CCCTTGCTTTTGTGT--------AGC +i HLfalTin1.ML997056 I 14 C 0 +s falPer1.KB390927 1430798 18 - 2055886 -------CCCTTGCTTTTGTGT--------AGC +i falPer1.KB390927 I 14 C 0 +s HLstrHab1.CM013776 2519986 7 - 7441144 --------------------------AGATGGG +i HLstrHab1.CM013776 C 0 C 0 +s HLtytAlb2.CABFVR010005803 1318274 18 + 2295574 --------CAATTCTGTTCTAAGATG------- +i HLtytAlb2.CABFVR010005803 C 0 C 0 +s halLeu1.KL869536 4538528 25 - 5862070 --------CCTTTCTATTCTAAGGTGGAAAAGA +i halLeu1.KL869536 C 0 C 0 +s aptFor1.KL225800 1006405 21 + 2472740 ------------TCTGTTCTAAGGTGGAAAAGA +i aptFor1.KL225800 C 0 C 0 +s HLcalAnn5.CM012136 503907 25 + 6461084 --------CCTTGCTATTCTAATGCACAAAAGA +i HLcalAnn5.CM012136 C 0 C 0 +s cucCan1.KL448280 2754948 25 - 4414095 --------tctttcttttctatggTGGAAATGA +i cucCan1.KL448280 C 0 C 0 +e HLcalPug1.KQ482592 384241 439 + 1956262 I +e HLlorGal1.manual_scaffold_27 1988627 80 + 7065481 I + +a score=231678.000000 +s galGal6.chr6 20883647 118 + 36374701 cagcactgggaaggaaaatgatttctgatAAAGACTTGGAGCCACAGAACGCTT-GGactagttttgcttttgcaaacagcTTAGTTTCAGGCTAAATTGAACCTGatcactgctttaa--------------------------------------------------------------------------------------------------------------------------------- +s falChe1.KB397934 1119966 118 + 4640569 CAGCATCGGGAAGGAACGTGCTTTCTGAGAAAGCCTCAGGTCTGTAGCGCACTT-AAACCAGTTTTGCTTGTGCACGCGGCTTAGTTTCAGACTAAATCAGACCTGACAACTGCTGTAG--------------------------------------------------------------------------------------------------------------------------------- +i falChe1.KB397934 C 0 C 0 +s HLfalTin1.ML997056 33510701 118 + 37096019 CAGCATCGGGAAGGAACGTGCTTTCTGAGAAAGCCTCAGGTCTGTAGCGCACTT-AAACCAGTTTTGCTTGTGCACGTGGCTTAGTTTCAGACTAAATCAGACCTGGCAACTGCTGTAG--------------------------------------------------------------------------------------------------------------------------------- +i HLfalTin1.ML997056 C 0 C 0 +s falPer1.KB390927 1430816 118 - 2055886 CAGCATCGGGAAGGAACGTGCTTTCTGAGAAAGCCTCAGGTCTGTAGCGCACTT-AAACCAGTTTTGCTTGTGCACGCGGCTTAGTTTCAGACTAAATCAGACCTGACAACTGCTGTAG--------------------------------------------------------------------------------------------------------------------------------- +i falPer1.KB390927 C 0 C 0 +s HLstrHab1.CM013776 2519993 99 - 7441144 TATCACAGGGAG-------------------AGCCTTAGTTCTACAGCACACTT-AAACTAGTTTTGCTTGTGCACACAGCTTAGTTCCAGACTAAATCAGACCAGACTATTTGAGCAA--------------------------------------------------------------------------------------------------------------------------------- +i HLstrHab1.CM013776 C 0 C 0 +s HLtytAlb2.CABFVR010005803 1318292 111 + 2295574 -------GGGAAGGAAAATGCTTCCCGAGAAACCCTTATGTCTCCAGCTCACTT-GAGCTAGTTTTGCTTTTGCACTCAGCTTAGTTTCAGACGAAATCAGACGTGACAACCGCTGTAA--------------------------------------------------------------------------------------------------------------------------------- +i HLtytAlb2.CABFVR010005803 C 0 C 0 +s halLeu1.KL869536 4538553 107 - 5862070 CAGCATTGGGAAGGAAAATGCTTTCTCAGAAAGCCTGAAGTGTAAAACACATTT-AAACTAGTTTTGCTTTTGCACACAGCGTAG-----------ATCAAACCTGACAACTGCTGTAA--------------------------------------------------------------------------------------------------------------------------------- +i halLeu1.KL869536 C 0 C 0 +s aptFor1.KL225800 1006426 119 + 2472740 CGGCATCGGGAAGGAAAATGCTTTCTGAGAGCGCCTTAGGTCTACAGCACACTTAAAACTAGTTTTGCTTTTGCAAACAGCTTCGTTTCAGACTAAATTGGACCTGACAACTGCTGTAA--------------------------------------------------------------------------------------------------------------------------------- +i aptFor1.KL225800 C 0 C 0 +s HLcalPug1.KQ482592 384680 111 + 1956262 cagtgttgggaaggaaaatgctttctgagaaagccttaggtctagagcacactt-aagctagttctacttttg-------gttggtttcagactcgatcggacctgacaattgctgtaa--------------------------------------------------------------------------------------------------------------------------------- +i HLcalPug1.KQ482592 I 439 I 708 +s HLlorGal1.manual_scaffold_27 1988707 49 + 7065481 --------------------------------AATAAAGAGCTACTGC-CAACT-ACATTGGTTTTGTGTGTA---AGAGTTGAGC------------------------------------------------------------------------------------------------------------------------------------------------------------------ +i HLlorGal1.manual_scaffold_27 I 80 C 0 +s HLcalAnn5.CM012136 503932 198 + 6461084 TGGCACTGGGAGGGGAAATCCTTTCTGAAAAATCCTTAGGTCTACAGCACACTT-AAACTACCTCTGCTTTTGCACCGACTTTAGTTTCAGACTAAATGGGGCCTCACAACTGCTGTGAAAGAGTTAGAGCAAAGATAATTAACAGTTCCTCTAACTGCAGGATTTTTGTGTTCAGAGT---------------------------------TCACCAG--------------TGAAG--TAGTTTAA +i HLcalAnn5.CM012136 C 0 C 0 +s cucCan1.KL448280 2754973 214 - 4414095 CAGCATCgggaaggaaaatgctttcttagAAAGCCTTAGGTCGACAGCACACTT-AAACTACTTAT-CCTTTGCAAACAGCTTCGTATCAGACTAAATCAGACCTGGCAACTACTCTAAAGGAGTGTGAGCAAACTTAATTAAC--------------------------------AGTTAGTGTTAACTacattaattttgtgtttagaGGTCACCAGGTCAGCACACTGCCTTAAGGCTGCTTTAT +i cucCan1.KL448280 C 0 C 0 + +a score=-20483.000000 +s galGal6.chr6 20883765 10 + 36374701 aa----------------------------------------------a-gag----------------tTTG +s HLlorGal1.manual_scaffold_27 1988756 5 + 7065481 ------------------------------------------------A-GG-------------------TG +i HLlorGal1.manual_scaffold_27 C 0 I 48 +s HLstrHab1.CM013776 2520092 52 - 7441144 -ACACAATAAAGAGCTAACTACATTGGTTTTGTGTGTAAGAGTTAACCA-GG-------------------TG +i HLstrHab1.CM013776 C 0 I 59 +s falChe1.KB397934 1120084 9 + 4640569 -A----------------------------------------------A-GAG----------------TTTG +i falChe1.KB397934 C 0 I 120 +s HLfalTin1.ML997056 33510819 9 + 37096019 -A----------------------------------------------A-GAG----------------TTTG +i HLfalTin1.ML997056 C 0 I 121 +s falPer1.KB390927 1430934 9 - 2055886 -A----------------------------------------------A-GAG----------------TTTG +i falPer1.KB390927 C 0 I 120 +s HLtytAlb2.CABFVR010005803 1318403 10 + 2295574 -A----------------------------------------------AGGAG----------------TTTG +i HLtytAlb2.CABFVR010005803 C 0 I 121 +s halLeu1.KL869536 4538660 9 - 5862070 -A----------------------------------------------A-GAG----------------TTTG +i halLeu1.KL869536 C 0 I 121 +s aptFor1.KL225800 1006545 9 + 2472740 -A----------------------------------------------A-GAG----------------TTTG +i aptFor1.KL225800 C 0 I 121 +s HLcalAnn5.CM012136 504130 19 + 6461084 AA----------------------------------------------A-GAG-------CTGTCACCTCCTG +i HLcalAnn5.CM012136 C 0 C 0 +s cucCan1.KL448280 2755187 26 - 4414095 AA----------------------------------------------A-CAGCACCAAACCCTCACCTCCTG +i cucCan1.KL448280 C 0 C 0 +e HLcalPug1.KQ482592 384791 708 + 1956262 I + +a score=2587130.000000 +s galGal6.chr6 20883775 498 + 36374701 GACAGTGAAAGTTTGACACTCACCCCTGTCCATCT-GTTCATTCCTAGGTGAATCTAAGCAGTTTTGATCCAGAAGACCCGGGACCAGCAGCGATGCACTGTGTGCTTCCTTCGATGGAGAGCTCTTAGGAAGCATGAATCGCTGCAGTTCTGATGACAGCCTTCCACATG---ACAACAACAATGGGGATAAGGAGCTCATCCACGTTCCTGCCCTCATTTCTCCTGGCTCCGCTGATGATGTTGACCTGAGCCCGCCAGACATCACAGTCGCTAGCCTGGATTGTGACCCAATGCCTTTCCAGTGCAGCCCACCCCAAGCAGAGCCCAAGTGCCTGAACAGCAGCACCTCCTTGA-GGAGCAGGTCGGCATCGCTGAGGAGAAGCCGAGCTTCTTGGAGGCGGACTTAGAGCCAGGCAGCCAGTCCCAggcaccaggcagcagcacaagctCTGAACTAGTCTCTCAATTCCCAGAG-----GAGTCCCTCTTTTACTCTGGATCT +s HLlorGal1.manual_scaffold_27 1988809 503 + 7065481 GACAATAAGC-TTTGACACAAACCCATCTCTGTCT-ACCCTTTTCTAGGTGAATCTAAACTCTTTAGGCCCAGGAGACCTAGGTCCAGTAGCGATGCATTGTCTGCTTCCTTCAACGGAGAGCTCTTGGGAAACATGAATCGCTGCAATTCTTATGACAACCTTCCTCACG---ATGACAGCGATGGGGATGAGGGGCTCATCCACGTTCCTGCCCTCATTTCTCCTCGCTCTGCTGAAGATGTTGACCTGAGCCCACCCGATATCGGAGTCGCAAGCCTGGATTTTGACCCGATGTCTTTTCAGTGCAGCCCACCCCAAGCAGAGTCCGAGTGCCTGGATAGCAGCACCTCCCTGCTGGAGTCCATTGGCATAAATAAGGAGAAGCCAAGCCTCCTAAAGAAAGATGTAGAATCAGGCAGCCAGTCGCAGACACCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCACGAGAAGATGAGTCCCTTCTTCACTCTGGACCT +i HLlorGal1.manual_scaffold_27 I 48 N 0 +s HLstrHab1.CM013776 2520203 503 - 7441144 GACAATAAGA-TTTGACACTAACCCCTCTCTGTCT-ACCCTTTTCTAGGTGAATCTAAACTCTTTAGGCCCAGAAGACCTAGGTCCAGTAGTGATGCATTGTCTGCTTCCTTCAACGGAGAGCTCTTGGGAAACATGAATCGCTGCAATTCTTATGACAACCTTCCCCACG---ATGACAGTGATGGGGATGAGGGGCTCATCCACGTTCCTGCCTTGATTTCTCCTCGATCTGCTGAAGATGTTGACCTGAGCCCACCAGATATCGGAGTCGCAAGCCTGGATTTTGACCCTATGTCTTTTCAGTGCAGCCCACCCCAAGCAGAGTCCGAGTGCCTGGACAGCAGCACCTCCTTGCTGGAGTCAGTTGGCATAAACAAGGAGAAGCCAAGCCTACTAAAGAAAGATTTAGAATCAGGCAGCCAGTCCCAGACACCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCATGAGAAGATGAGTCCCTTCTTCACTCTGGACCT +i HLstrHab1.CM013776 I 59 I 15517 +s falChe1.KB397934 1120213 507 + 4640569 GACAGGAAGATTTTGACACTAACCCACGTGTATCT-ACTCGTTCCTAGGGGAGTCTAAACTGTTTCGACCCAGAAGACCTAGGTCCAGTAGTGATGCGTTGTCTGCTTCCTTCAACGGAGAGCTCTTAGGAAACATGAACCGCTGCAATTCTTACGACAACCTTCCCCACGACAATGACAGTGATGGGGACGAGGGGCTCATTCACGTTCCTGCCCTGATTTCTCCTCGGTCTGCTGAAGATGTTGACCTGAGCCCACCGGATATCGGAGTCGCCAGCCTGGATTTTGACCCAATGTCTTTCCAGTGCAGCCCTCCCCAAGCAGAGTCCGAGTGCCTGGACAGCAGTACCTCCCTGCTGGAGTCAGTTGGCATAAATAAGGAGAAGCCAAGCCTGCTAAAGAAGGATTTAGAATCAGGCAGCCAGTCCCAGACCCCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCAAGAGAAGATGAGTCCCTTCTTTACTCTGGACCT +i falChe1.KB397934 I 120 N 0 +s HLfalTin1.ML997056 33510949 507 + 37096019 GACAGGAAGATTTTGACACTAACCCGTGTGTATCT-ACTCGTTCCTAGGGGAGTCTAAACTGTTTCGACCCAGAAGACCTAGGTCCAGTAGTGATGCGTTGTCTGCTTCCTTCAACGGAGAGCTCTTAGGAAACATGAACCGCTGCAATTCTTACGACAACCTTCCCCATGACAATGACAGTGATGGGGACGAGGGGCTCATTCACGTTCCTGCCCTGATTTCTCCTCGGTCTGCTGAAGATGTTGACCTGAGCCCACCGGATATCGGAGTCGCCAGCCTGGATTTTGACCCAATGTCTTTCCAGTGCAGCCCGCCCCAAGCAGAGTCCGAGTGCCTGGACAGCAGTACCTCCCTGCTGGAGTCAGTTGGCATAAATAAGGAGAAGCCAAGCCTGCTAAAGAAGGATTTAGAATCAGGCAGCCAGTCCCAGACCCCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCAAGAGAAGATGAGTCCCTTCTTTACTCTGGACCT +i HLfalTin1.ML997056 I 121 N 0 +s falPer1.KB390927 1431063 507 - 2055886 GACAGGAAGATTTTGACACTAACCCACGTGTATCT-ACTCGTTCCTAGGGGAGTCTAAACTGTTTCGACCCAGAAGACCTAGGTCCAGTAGTGATGCGTTGTCTGCTTCCTTCAACGGAGAGCTCTTAGGAAACATGAACCGCTGCAATTCTTATGACAACCTTCCCCACGACAATGACAGTGATGGGGACGAGGGGCTCATTCACGTTCCTGCCCTGATTTCTCCTCGGTCTGCTGAAGATGTTGACCTGAGCCCACCAGATATCGGAGTCGCCAGCCTGGATTTTGACCCAATGTCTTTCCAGTGCAGCCCTCCCCAAGCAGAGTCCGAGTGCCTGGACAGCAGTACCTCCCTGCTGGAGTCAGTTGGCATAAATAAGGAGAAGCCAAGCCTGCTAAAGAAGGATTTAGAATCAGGCAGCCAGTCCCAGACCCCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCAAGAGAAGATGAGTCCCTTCTTTACTCTGGACCT +i falPer1.KB390927 I 120 N 0 +s HLtytAlb2.CABFVR010005803 1318534 504 + 2295574 GACAGACAGATTTTAACACTAACCCGTGCTTCTCT-GCTCCCTCCTAGGTGAATCTAAGCTGTTCCGACCCAGAAGACCCAGGTCCAGTAGCGATGCGTTGTCTGCTTCCTTCAACGGAGAGCTCTTGGGAAACATGAACCGCTGCAATTCTTACGACAACCTTCCCCACGAC---GACAGCGACGGGGACGAGGGGCTCATCCACGTTCCTGCCCTGATTTCTCCTCGATCTGCTGAAGACGTTGACCTGAGCCCACCGGACATCGGAGTCGCTAGCCTGGATTTTGACCCTATGTCTTTTCAGTGCAGCCCGCCCCAAGCAGAGTCCGAGTGCCTGGATAGCAGCACCTCCCTGCTGGAGTCAGTCGGCATCAACAAGGAGAAGCCAAGCCTACTAAAGAAGGATTTGGAATCGGGCAGCCAGTCCCAGACACCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTACCAAGAGAAGATGAGCCCCTTCTTTACTATGGACCT +i HLtytAlb2.CABFVR010005803 I 121 N 0 +s halLeu1.KL869536 4538790 507 - 5862070 GACAGTAAGATTTTGACACTAACCCGTACGTATCT-ACTCATTCCTAGGTGAATCTAAACTGTTTCGACCCAGAAGACCTAGGTCCAGCAGTGATGCATTGTCTGCTTCCTTCAATGGAGAGCTCTTAGGAAACATGAATCGCTGCAATTCTTACGACAACCTTCCCCACGACCATGACAGCGATGGGGATGAGGGGCTCATCCACGTTCCTGCCCTGATTTCTCCTCGATCTGCTGAAGACGTTGACCTGAGCCCACCGGATATCGGAGTCGCTAGCCTGGATTTTGACCCAATGTCTTTTCAGTGCAGCCCACCCCAAGCAGAGTCCGAGTGCCTGGATAGCAGTACCTCCCTGCTGGAGTCAGTCGGCATAAATAAGGAGAAGCCAAGCCTAGTAAAGAAGGATTTAGAATCAGGCAGCCAGTCCCAGACACCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCAAGAGAAGATGAGCCCCTTCTTTACTCTGGACCT +i halLeu1.KL869536 I 121 N 0 +s aptFor1.KL225800 1006675 507 + 2472740 GACAGTAAGGTTTTGACGCTAATCCGTATGTATCT-ACTCGTTCCTAGGTGAATCTAAACTGTTTCGACCCAGAAGACCTAGGTCCAGTAGCGACGCATTGTCTGCTTCCTTCAATGGAGAGCTCTTAGGAAACATGAATCGCTGCAATTCTTACGACAACCTTCCCCACGACAACGACAGTGATGGGGACGAGGGGCTCATCCACGTTCCTGCCCTGATTTCTCCTCGATCTGCTGAGGATGTGGACCTGAGCCCACCGGATATCGGAGTCGCTAGCCTGGATTTTGACCCAATGTCTTTTCAGTGCAGCCCGCCCCAAGCAGAGTCCGAGTGCCTGGATAGCAGTACCTCCCTGCTGGAGTCAGTTGGCATAAATAAGGAGAAGCCAAGCCTACTAAAGAAGGATTTAGAATCAGGCAGCCAGTCCCAGACACCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCAAGAGAAGATGAGTCCCTTCTTTACTCTGGACCT +i aptFor1.KL225800 I 121 N 0 +s HLcalPug1.KQ482592 385499 505 + 1956262 gagagcaagatcttgACCCTAATCCACATGTCTCTAACTCGTTCCTAGGTGAATCCAAGCTGTTCCGACCCAGAAGACCCAGGTCCAGCAGCGATGCCCTGTCTGCATCCTTCAACGGAGAGCTCTTAGGAAACATGAACCGCTGCAACTCCTACGACAACCTTCCCCACG---ACGACAGCGACGGGGACGAGGGGCTCATCCACGTTCCTGCCCTCATCTCTCCTCGATCCGCCGAAGACGTCGACCTGAGCCCCCCCGACATCGGCGTCGCCAGCCTGGACTTTGACCCCATGTCTTTCCAGTGCAGCCCACCCCAGGCGGAGTCCGAGTGCTTGGAGAGCAGCACCTCCCTGCTGGAGTCCTCTGGTCTGAACCAGGAGAAGCCCAGCCTTCTAAAGAAGGATTTGGAATCAGGCAGCCAGTCCCAGACACCGGGCAGTGCCACGAGCTCTGAGCCGGTCTCCCCTTTCCAAGAGAAGATGAGTCCCTTCTTTACGCTGGACCT +i HLcalPug1.KQ482592 I 708 N 0 +s HLcalAnn5.CM012136 504149 503 + 6461084 CACAGTAAGCTGTTGGCACCAGCCCATATTT-TCT-ACTTGTTCCCAGGTGAATCGAAGCTGTTCCGACCCAGAAGGCCTCGGTCCAGCAGCGACGCGCTGTCCGCTTCCTTCAACGGAGAGCTCCTGGGGAACATGAACCGCTGCAATTCTTATGACAACCTTCCCCACG---ATGACAGTGATGGGGACGAGGGGCTCATCCACGTTCCTGCCCTCATCTCTCCTCGATCTGCTGAAGATGTTGACCTGAGCCCCCCGGATATCGGAGTTGCCAGCCTGGATTTTGACCCCATGTCTTTTCAGTGCAGTCCACCCCAAGCGGAATCCGAGTGCCTGGACAGCAGCACCTCCCTGCTAGAGTCAATCGGAATCAATAAGGAGAAGCCAAGCCTCCTAAAGAGGGATTTTGAGTCAGGCAGCCAGTCCCAGACACCAGGCAGTGCCACGAGCTCTGAGCCGGTCTCCCCTTTCCAAGAGAAGATGAGTCCCTTCTTTACTCTGGACCT +i HLcalAnn5.CM012136 C 0 N 0 +s cucCan1.KL448280 2755213 499 - 4414095 GACAGTAAG--ATTGACACTAACTCATATGTATCT-ACTCATTCCCAGGTGAATCCAAACTGTTTCGACCCAGACGGCCTCGGTCCAGCAGTGACGCATTGTCTGCTTCATTCAACGGAGAGCTCTTAGGAAACATGAATCGCTGCAATTCCTACGACAACCTTCCCCACG---ATGAAAGTGATGGGGACGAGGGGCTCATCCATGTTCCTGCCCTGATTTCTCCTCGATCTGCTGAAGACGTTGACCTGAGCCCACCAGATATCGGAGTCGCAAGCCTGGATTTCGACCCAATGTCTTTTCAGTGCAGCCCACCCCAAGCAGAGTCCGAGTGCTTGGA---CAGCACCTCCCTGCTGGAGGCAGTCGGCATAAATAAGGAGAAGCCAAGCCTGCTAAAGAAGGATGCAGAATCAGGCAGCCAGTCCCAGACACCAGGCAGTGCCACGAGCTCTGAGCCAGTCTCCCCTTTCCAAGAGAAGATGAGTCCTTTCTTTACTCTGGACCT +i cucCan1.KL448280 C 0 N 0 + +a score=0.000000 +s galGal6.chr6 20884273 485 + 36374701 CTCTTGCCCTTCTGTCTTTTAATGAGTAAAGGCATACTCATGTTCAATCAACACCTTGACTCCGTTCTGGGTCTTCATTCACAGGGACTGGGTTCAATATAAAAAACTGCAAGGCACTTCAGCTTGCAATGGTATGGCttttggctgtgtgctgcagagctctggcttTGTTCCCACGTTGTATGGCTTtggggtgtctgcagggagagaTCTGCTACTGTTTTGTATGGTAGCAAATAGCATCCTGGTTGCTGCTCATTAAAATCCCCAGCTTCCTGGCTGCAGCTTTGTGACGCTTTTCAATGCTTTGTGATGCTCTTTGATGCCTTGGCTGTTCCCTgacacctctctgctcctgcagtgcccaTGCCTGGCTCTCACCTATGAGCTGGACTTGCAGATCATCTCCAGACCTGTCATCatcacactgcagcaggagaagtTCATCCTGACCAACGTCAGGCAGTACTGCTTGGTCTGCACGT +e HLstrHab1.CM013776 2520706 15517 - 7441144 I + +a score=7715.000000 +s galGal6.chr6 20884758 144 + 36374701 TACTGGCTGCCTCTGCATCTCTTCAGCTGCAGTACCATCTCTCAGGGAAGTCAGATGCAGTCACGCTCCTCTCTGAATTCTGTGCTTGATTAAATCCAAAGCAAATCTTTACACTACAATGAGCAAACCATGCTGTCTACCAAA +s HLcinPul1.manual_scaffold_3 87619818 135 - 117014874 TACTGACTGCCTCCCCATCTCTTTAAAAATAGTA---------AGGGAAGCCAGATGCAGTCACATTCCTTTCTGAATTTAGTGCTTGATTAAATCATAAGCAAATCATCAGACCATGATGAGCAAACCACCCTCTATACCAAA +i HLcinPul1.manual_scaffold_3 N 0 N 0 +e HLstrHab1.CM013776 2520706 15517 - 7441144 I + +a score=0.000000 +s galGal6.chr6 20884902 1155 + 36374701 GAGGCAGGTTTGTGCCCCAGGATTTCAGGAAGAAGGACTGGTATCAGCCCAGGGCAGGATCTGGgtcccctggggacagcagctaCAGCCCCTGTCCAATCCGCACGTGACCCTGAGCCATCCCTAGGGTGCCCTgtgacctcctcctcctcctctgcacagccccacaaaGCACAAGGGAGGTCTGGGAGGTGGCAATGCCTAGATGCCCCTGGAACAGGCCAGGAATGGTgccctccatcccatccatccccaggcCCTGATGCCCAGTCTCAGGGCTCCAGGATCTCTCAGTGtccctgttccagggcctcaccacccttactacagaatcacagaatggtgcTTATCatagaaaacttcttccttatgaaGAATCTAAAGAGTCCGTTCCCTCCTTTGTTGCAGCCCCCCTTGAGATACtgaaggccactctcaggtcttcctggagccttctcttctccaggctgcacagcctcagctctcttaGCCGTCCTTGGAAGGTGGTATTCCATCCCTGggaccatttttgtggccctcttctggtcACGccccaacaggtctatgtctctcctgtactgaggactccacatctggatgcagtgctcacagtgcagaggggcagatcccctccctcaccctgctggccgCGCTGCTTCGGATGCTGCCCAGGGTATGGCAGGCTTTCTGGGCTGagggcactctgctggctcatgcccagttgccatccatcagtaccctcaggtcctttttggcagggctgtgctccatccttacatcccccagcttgtactgaaaGCAGGGGgctgccatgacccaggtgcaggccTTGCACCCGGATTCATTCAACCTCATGAGCTTCctctgggcccactgctcagcatgTCTCTCTGTTCGCATCCCATCCTTTGGACGTGCGGGCTGCACCATGCAGCTTGGGGTCACCTGCAAATCTGAGGGGGCACTCAATCCCAGTGTTggatgtcattgatgaagatgttaaagagcactggttccAGTACTGGACcacctgagggacaccactggGTCTCCACCCGGACATGGAGCCCTTGACCACCACTCTGTGCATGCGATCTTGCAGCCAGTCCATCCAACAgtctacccatcaaatccatatcttccCAATTTGGAGCT +e HLstrHab1.CM013776 2520706 15517 - 7441144 I + +a score=9865.000000 +s galGal6.chr6 20886057 303 + 36374701 GCTCTTTTTGGTTCCGATCTGCGGGGCCTGcaaatccacacagtccaataa--ctGACGATCCCTTTCCTCCACCGGGCTGGAGGCAGGGCCCCTCTAGT-CTGCCTATTGTATTCATCGCTCACTTCTCATGAACGTGAACCAGTACACTTAACCCCAGAAATCAGAT------CCGACCATCCAGTCTGTCCGGAGAGCTGCCCATTAGAAACTggagctgtcatttttttttctggaagaatccccttttcctattgttttctctgcaactcACCTGCCCGTGCATGTAGGACTGAGGTGGGTTTTCCA +s HLstrHab1.CM013776 2536223 305 - 7441144 GCTGCCTTTGACTCCG---TGCCAGCCCATCTAATCCACATAGTCCAGCAAACCTGGTTGCCCCTCTCCTGCACGTGGCTGGAGGCAGGGCCACTCTAGTCCTGGTTGTAGGATGCTCTGCTTCCTTCTTGTAAAAGTGGGTTAGAATTCATCACCACAG-GGTCAGATGTAAGGTCACCCTCCCTCTCTGCATGGGAAGCTGCCCACTGGAGACTGGAGCAGGGACTCTTCTGT---GCAAACCCCCATTTGTGGATGTTTTTCCTTGCAATTCACATACCCGTGTGTCCAGGACTGAGGTAGGTTTTCCA +i HLstrHab1.CM013776 I 15517 N 0 + +a score=116027.000000 +s galGal6.chr6 20889624 686 + 36374701 ACCTCCAAACTCATCtaggagaataatttttttcgtaatattcaacctgagagtgctgaggtatAATGC--TCCACAGTCAGACCTACATGTGACTGTGGGGCCCTGATGTTCTTCCCATGGGCACACAAA--CACTCAATGaccatactgggaggcaaaaagcacagcagaagacGATCCTCAGCTGCAGACGCTGATGACA------------------------------------------------------CTCATAGCAATAAGAGAGAAGTATGTGTGAAAGGCTTACGTTACAAAAGGAATGACTGATTGTCAGCAGAGATGTGGTCTAGCACTTGTTCTG-------CTGGAGCAACTCCAGTGTGGCCATAGATTTGCACTACGGAGCAAAAACAAATCTCAGCCGGAATAATACTGTCAAGCCTAAA-ACCAACATTCGGCTCCTAAAATCCTGAGTTGTtggttgaaaaatatttcagcccTGATGGTGTGCTCAGTGGAGTTGTTTCTGTGCATATTGTCTCC-----------GTTCCCTCTActgcaaacaaaatgcattcaaatcCCCAGGAAGAAGCTCATGAATACTGCCGCACCCTTCCCGATCAATCCTTCCAGACAGCTCGGCTTTTTCTAAGC-AAAACATACCTATGGAAAAGTCTCAGG-GGGTAACCTAGAAGCCAGCGAACGCTGACCTCTTCTGGATACTTCTAGAAATGGCATCGAGGCATAA-----CCATTGCAGCATA-------------------------------------------------------------------------- +s halLeu1.KL869074 15614759 768 + 16075551 ACATCCTAACTAGGCACTGAGAAGTACTTTTTTTA-AATGT---ACCTGGGAGTGCTGAGGTATAACGCAGGCGAGTGTCGCACCCATGTACAACTATGGGGCCCTGATGTTCTGGCCACAGGCACACAAACTCTGTCAGTGACCATTTTGGGGAGC-AAAAGCACAGCAGAAGCCGATCCTTGGTTGCAGACGCTGATGATATTCTTTATTCCAGAAGGCACGAGGCTCGGCTGCACTGCTGCAATCCCAGTGACCTTCACAGGGATCAGAGAGAACTGAGCCTGAAAGCTTTAAGTGAAAAAGCTAATCAGTGATTGCCAGCAGAGACATTGTCTAGTACTTCTCCTGTGTTGGCCTGGAGCAACTCCAGGGCGCTCACAGATTTGTCCCACT--GCAGAGAGAAACTTCAGTCTGACTAATACTGTCAAGCCTAAACACCAGATTTTACCTGTTAAAATCCCATGACGCTAGGTTTAAAATATTCAA----TGATGATGTGGTCAGCAGAGCTGTTTGCCTGCATATTGTCTCCCGGATCCCTGGTTCCCCTATACTGCTAAGAAAATGCATTCAAATCCCCAGGAAGAAGCACATTAATACAGCTGCATCCCTCTTGATCAAACTTTCCAGACAGCTCAGCATTTTCTAAATAAAAACATACGTTTGGAAAGGTTTCAGGCAAGTAACCTAAAAGCGAGGGCACACTGACCTCTTCTGGACATTTCCAGAAGTGGCAACAAGGTAGAAAAGAGCTTGCACAGCTTTCTGTGCT----------------------------------------------------------TG------- +i halLeu1.KL869074 N 0 C 0 +s aptFor1.KL225833 568392 833 - 667431 ACATCCTAACTAGGCACTGAGAAGTACTTTTTTTA-AATGT---ACCTGGGAGTGCTGAGGTATAATGCAGACGAGCGTCAGACCCATGTACAACTATGGGGCCCTGATGTTCTGGCCACAGGCAGAAA----CTGTCGGTGACCGTTTTGGGGAGC-AAAAGCACAGCGGAAGCCGATCCTCGGTTGCAGATGCTGATGATATTCTTTATTCCAGAAGGCACGAGGCTCGGCTGCACTGCTGCAATCCCAGTGACCTTCATAGGGATCAGAGAGAACTGAGCCTGAGAGCTTTAAGTGAAAAAGCTAATCAGTGATTGCCAGCAGAGACATGGTCTAGTACTTCTCCTGTGTTAGCCTGGAGCAACTCCAGGGTGCTCACAGATTTGTCCCACT--GCAGAGAGAAACCTCAGTCTGACTAAGACCATCAATCCTAAACACCAAATTTTAGCTGTTAAAATCCCATGACGCTAGGTTTAAAATATTCAAGCCCTGATGATGTGGTCAGCAGAGCCGTTTGCCTGCATATTGTCTCCCTGATCCCTGGTTCCCCTGTACTGCTAAGAAAACGCATTCAAATCCCCAGGAAGAAGCACATTAATACAGCTGCATCCCTCTTGATCAAACTTTCCAGACAGCTCAGCATTTTCTAAATAAAAACATACGTTTGGAAAGGTTGCATGCGATTAACCTAAAAGCGAGGGCACACTGACCTCTCCTGGGCATCTCGAGAAGTGGCAACAAGGTAGAAAAGAGCCTGCACAGCTTTCTGTGCTTGGTTATACATAACTATCCAGGGGCAACAAAGCAATGGTCGTATTTTCCTGTATTTAGTGCATAATA +i aptFor1.KL225833 N 0 I 165 + +a score=-7300.000000 +s galGal6.chr6 20890310 23 + 36374701 --------------------------------------------------CTTCTACGCAA------CGAGGCGTAACC--------------------------------------- +s halLeu1.KL869074 15615527 118 + 16075551 GTTATACGTAACTATCCAGGGGCGACAAAGTAAAGGTCATATTTTCCTGTATTCAGTGCGTAATAGTTAAGTCTTCTCCTCCAGTTATTTTTATCTTGATAAGGATGTTCCCTTAACC +i halLeu1.KL869074 C 0 I 108 +s aptFor1.KL225833 569390 53 - 667431 ------------------------------------------------CTTTTAAGCGCAT--------AGTCGTGGCCCACAAACACTACTAGTATGACCAGAA---------AAAC +i aptFor1.KL225833 I 165 C 0 + +a score=81546.000000 +s galGal6.chr6 20890333 537 + 36374701 GATCCCTTTT-----------------------------------------------------ATTCACCTCCCTCAAGAAAAGGTCAACTCTCAGCCCCACTTTCATCTGCTTGGTGAGCTGGCTCAGTTCAGAGTTCTGCAGGAGAGCAAACAAGGTGCAGAAAGAAAGACTCCAGAGGTTACAGCTATGTTTCAGCTG-TTAAGACTTACATCTTTTC----CTTCTATGCTC--ATTTGCTCTTGT------------------------------------------------------------------------------------GCCTATGTG----GAGTAGAAACAGATTGAAATAATTGAACTGATTGAAATCTTAAGGGTCATGCATGTTTGATGCAATACTGTGCAGATCCTAAATTGATACGGCTTATGGCAAGTAATGTAATTTACCACCAGGTCAAAAAAGAATTATGTTTCTTGAAACATATGATTTTGCACTACGGATATACACTATTGACTAGAGAAGGTATCATGACATCAAGACAAAGCTGCAAATGAAGCTCTCACTGAAAAGGCATTTATAAGCTGTGTTTCAGATGTCAGGGTCTTGAGTG--CCCTCCCACCTTTGTT----------------------CCTGCAGCTAGCA-CTTTTCCTCACCAGCTACCACCTCCACCT +s halLeu1.KL869074 15615753 690 + 16075551 GATTTCTTTTAAGTGCATAGCAGTGGCCCACAAACACTACTAGCATGACCAGAAAAGCCATCACTTTACTACCCCCAAGAAAAGGTCAACTCTCAAGCATACTATCATGTAGCTAGTAATTTGCCTAGGTTCAGAGTTCTGCAGAGGAGCAAA---------------AAGACTCCAGAACTCCTATCTGTATTTTACCTGTTTAAAAACTAAAAGTTTCCACCTCTTATGTGCTCCTGTTTGCCCTTGTAAGACATACTTACTAGCTCCAGCAATAGCAATCCTTTCTACAGCTAGAAGCCAACATCTTCCATCAATAGTAATGGCACTCAGGGCCTATATATATATAGTAGGCAGACATGGAAATAGTTGAGCTGATTGAGATCTGTAGGCTCATGCATGTCTGATGCAATACTATGCATAGCCTGGATTGACACTGCGTACATC----AGTAGGATTGACCACCAAGTC-AAAAAGAGCCATGTTTCTCAAGGAGCATGACTTTGCACTACAAATACACGCTGTTGTTTAAAGCAGGTATTGTTACATTAAGAGAAATTTGCAGATTGGGTTCTCACTGAGAAGGGATTTATCAGCTGTGTTTCAAGTGCCGGGGTCCTGAGTGTAGCCTCAGGCTTTCCTTGCCCCAAGAAGTCTCACCTGGACCTGTAGCCAGCCCCCTTTGCTCATCGACTGTGACTTCCATTT +i halLeu1.KL869074 I 108 I 9 +s aptFor1.KL225833 569443 633 - 667431 CATCACTTTT-----------------------------------------------------CTTTACCACCCCCAAGAAAAGGTCAACTCTCAAGCATACTATCACCTAGCTAGCAATTTGGCTACGTTCAGAGTTCTGCAGAGAAGCTAA------------AGGAAGACTCTAAAACTCATATCTATATTTTACGTGTTTAAAAACTAGAAGTTTCCACCTCTTATGTGCTCCTGTTTGCCCTTGTAAGACATACTCACTAGCTCCGGCAATGGCTACACTTCCTGCAGCTAGAAGCCAACATCTTCTATCAGTGGTAATGGCACTGATGGCCTATATA-------TAGACACACATGGAATTAATTGAGCTGATTGAGATCTGTAGGCTCATGCATGTTTGATGCAATACTACGCATAGCGTGGATTGATACTGCGTACATC----AATAGGATTGACCACCAAGTC-AAAAAGAGCCATGTGTCTCAAGGAGCATGACTTTGCACTACAAATACACGCTATTGTTTAAAGCGGGTATCGTTCCGTTAAGAGAAATTTGCAGATTGGGTTCTCACTGACAAGGGATTTATCAGCTGTGTTTCAAGTGTCAGGGTCCTGAGTGTACCCTCAGGCTTTCATTGCCCCAAGCAGTCTCACCTGGGCCTCCACCCAGCACCCTTTGCGCATCAACCGTGACCTCCATTT +i aptFor1.KL225833 C 0 I 319 + +a score=-4525.000000 +s galGal6.chr6 20890870 19 + 36374701 -------------CCCACTCG--------------------------------TGTTGCTGCCC---------------------------- +s halLeu1.KL869074 15616452 92 + 16075551 CCCAGGAGACTCGGCCGTTCTTTGAGGCAATTGTAGGATAGTCCTGGCTGTAGTTCTGTCTTTTGACTTTGCTGGGTAAGCTCTGGACTTGT +i halLeu1.KL869074 I 9 I 230 +s aptFor1.KL225833 570395 19 - 667431 -------------CCTGTTCT--------------------------------TGTTGCTTTTT---------------------------- +i aptFor1.KL225833 I 319 C 0 + +a score=136605.000000 +s galGal6.chr6 20890889 858 + 36374701 AGCACCAGGTAGCTCAT------------------TTTCCCCTGCTCCGA------CATCACCCA-TCACTTTCTGTGCTCCTTCCCTGAACAAATCAATG----TGTCATTTAGAAGTCCGCCTGAGCTGAGCAGTCacagctctgtttcctttccctgaAGCTG-GAGTACAGCCAGGCTACAGAGCAGCTAAGTAAGATAAAGTTGCATGGAAACAGATGTTaattaagaag-aagaaatgtaggCCACTCATCTCTTCTTCATAATGGTACGGGCTCTGATTGAGGTAAACATCCGTGAAATCCATCAGATGAGATTTACTTCTGCCAATGTGTTTAAGGCCGCTGGAAGCAATGGTGTAGCAGTCTCTGGCAtgcagctctggggaaggaggcCTGAGGAAATCCCACAGATAA--------TTCCACACTGGCCTAAGTAGCTTTCCC-ACCGGATGTctccctgtggggctgccaACACACTTCAGGCCTGGAGCTTAAACTGCATCAACTAGTCTTGA-----------------------------------------AAC----------------------------AG---------TGTCCTTAACGTAGCTAGGCTGTACTGGGAGAAAGAGTGTAATTTGGACTGAGCTGCATCCACACCAGGGGATGTCTTCTCAGCAAGGAATGGGAGCGGCCTACCAAAGGCAGGCAGTGCACGATGGTGGGAAGA---CGGTGTGTCTTCTGTCACAGCCTGTGCCTTCAGCTCTACTGGTGAtgtcccagcagagctgttttccCCAGTCCGGGTGGAATGGTGGGAAGTGATGCAACAGGAGCAGCGAGGATCAAGAGCAGGCgactgctgctgtgactgccAGGGGCTCATTGAAGGAGCACAGTTTCAGCCCGCTAAGGCCTCAGTCACCTTCACAACACC-CTTGAAAGGTGCTCGCTTACGCCGG +s halLeu1.KL869074 15616774 964 + 16075551 AGCACTAGATACATTATGTTAAGATGCTTAACATGTTCCCCCCACTCCAAGCACTACATCAGCCAGTGAGTTTCTATCATATTTCCATGAAAGAACAAATACCTCTGTGCCTTAGAAGTCCCCCTTAGCTGGGCAGTCATAGCTCTGTTTCCAATCC----ATCTGTAATTATAGCCAAATTACACAGCAGCTAAGTAAGATAAAGCTAGATAGAGGCAGATATCAATTAAGAAGAAAAAAATGAAAGCTACTCATTTCCTCTTCAGAACTATACGTATGCTAATTGAAGTAAACATTCATTAAGCCCATTAAATGAAATTTAGGTTTGATAATGCATTAAAGGCTGCTGGAAGCAATGGTATAGCAGTTCATAGCCTGCAGCTTTGGGGAAGGAAGGTTTAGGTTGTACCACAAATGATATGTAGCTTCCGCAAGGGCATTAGGAGTTCCCCTGACTGAGTATCAATCCATGGAGCTGGCAAATGCCCTAAGGCCTGGGGATTAAATTATATTAATTAGTCTTGATAGTGGGGGTTTGTTCAAGACAAAGTTATTATATTTGTACAAAC-------TGAATGCAATTATAGTGCTCTAGCCCATCTCTTGTCCTTAATGTAAATGGTTTATACCATGACAAGCATTGTAATTCTGACTGCCACGCATCCACACTAGGTGATGTCTTCCCAGCAAGTAATGGGAATGGCTTATCAAGGCCAGGCAAGCCTAGATCACGTGAAGACTCCTGTGTGTCTTCTGCCACAGCCTGCATCTTAAGCTCTGCTGGGGACGTTAAAGCAGACCTGTTTTCCCCAGTTCAGTTTAAACAGTGGAAGATGAGACAACA----CAACTAGGATCAAGAGCAAGCGACAACTTGTATGACTGCAAAGGGCTCTTTTGAAGGTAAAAGCGTCACCCCACAAGTACCT-AACCATCTTCACAACCCCACTATAGAGGTGCTTGCTGACTCCAG +i halLeu1.KL869074 I 230 I 3463 +s aptFor1.KL225833 570414 965 - 667431 AGCACTAGACACATCATGTTAAGATGCTTAACGTGCTTTCCCCACTTCAAGCACTGCATCAGCCAGTGAGTTTCTATCATATTTCCATGAAGGAATAAATATGCATGTGCCTTACAAGTCCCCTTTAGCTGGGCAGTCATAGCTCTGTTTCCATCCCCTAGAGCTG-GGTTATAGCCAAATTACAGAGCAGCTAAGTAAGATAAAGCTACAGAGAGGCAGATATTAATTAAGAAGAAAAAAATGAAAGCTACTCATTTCCTCTTCAGAACTATGTGTTTGCTAATTGAAGTAAACATTCATGAAGCCCATTAAATGAAATTTAGGTTTGATAATGCATTTAAGGCTGCTGGAAGCAATGGTGTAGCAATTCATAGCCTGCAGCTTTGGGGAAGGAGGGTTTGGGTTGTACCACAAATAATGTATAGCTTCCACAGGGGCATTAGGAGTTTCCCTGACTGAATGTCTCTCCATGGAGCTGGCAAATGCCCCAAGGCCTGAGGCTTAAACTTTATCAATTAGTCTTGA---CTTAGGTTTGTTTGAGGCAAAGTTATTATGTTTGTACAAACTATAATATGAATGCAATTATAGTGGTCTAGCCCATTTCTTGTTCTTAATGTAAAAAGGCTATACCATGACAAG-AGTGTAATTCTGATTGCACTGCATCCACACTGGGTGAGGTCTTCCCAGCAAATAAAGGGAATGGCTTATCAAGGCCAGGCAACCCTAGATCATGTGAAGACTCCTGCGTGTCTT------AAGCCTGTGCCTTAAGCTCTGCTGGGGATGTTAAAGCTGGCCTGTTTTCCCCAGTTCAGTTTGAACAGGGGAAGGTGAGACAACACAAGCAACTAGGATCAAGAGCAAGCGACAACTTGTATGACTGCAAGGGGCTCTTCTGGAGGGAAAAGTGGTACCCCACAAGTG----GGGTATCTTCAGAACCCCACTATAGAGGTGCTTGCTGACTCCAG +i aptFor1.KL225833 C 0 I 2965 + diff --git a/tests/real-excerpt.maf.block.idx b/tests/real-excerpt.maf.block.idx index 1fd879a..088e6a2 100644 --- a/tests/real-excerpt.maf.block.idx +++ b/tests/real-excerpt.maf.block.idx @@ -1,4 +1,4 @@ -# mafutils-index format=2 maf=real-excerpt.maf compression=none size=47241 mtime=1784644433.0 hash=md5:7146e0c45ac589e8998a3fe232ca5ad7 +# mafutils-index format=3 maf=real-excerpt.maf compression=none size=47241 mtime=1784644433.0 hash=md5:7146e0c45ac589e8998a3fe232ca5ad7 CM000994.3 125888838 82 84 135 13 17 1874 CM000994.3 125888920 187 223 275 9 1874 4416 CM000994.3 125889107 10 31 82 9 4416 5221 diff --git a/tests/real-excerpt.maf.scaffold.idx b/tests/real-excerpt.maf.scaffold.idx index dcdb6a9..a13c15f 100644 --- a/tests/real-excerpt.maf.scaffold.idx +++ b/tests/real-excerpt.maf.scaffold.idx @@ -1,2 +1,2 @@ -# mafutils-index format=2 maf=real-excerpt.maf compression=none size=47241 mtime=1784644433.0 hash=md5:7146e0c45ac589e8998a3fe232ca5ad7 +# mafutils-index format=3 maf=real-excerpt.maf compression=none size=47241 mtime=1784644433.0 hash=md5:7146e0c45ac589e8998a3fe232ca5ad7 CM000994.3 17 47241 diff --git a/tests/test_index_seq_lines.py b/tests/test_index_seq_lines.py new file mode 100644 index 0000000..a2ffecb --- /dev/null +++ b/tests/test_index_seq_lines.py @@ -0,0 +1,294 @@ +""" +Index format 3: the per-block sequence count (block.idx column 6) and the +seq_lines / max_seqs header totals count 's' lines only. Formats 1-2 also +counted i/e/q lines, which inflated the count on MAFs that have them. + +Fixture: multiz45way-excerpt.maf -- the `##maf` line plus 14 consecutive real +blocks (galGal6 chr6:20883553-20891747) of a 45-way multiz MAF, with 72 s, 58 i +and 12 e lines. It includes the block at chr6:20884758 (2 s + 1 i + 1 e lines) +that a format-2 index reported as 4 sequences. + +Malformed inputs for the error tests are built from copies of the fixture in +tmp_path; nothing here is invented sequence data. +""" + +import csv +import gzip +import os +import subprocess +import sys + +import pytest + +from mafutils.lib import bgzf +from mafutils.lib import common as COMMON + +TEST_DIR = os.path.dirname(os.path.abspath(__file__)) +REPO_ROOT = os.path.dirname(TEST_DIR) +MAF_FILE = os.path.join(TEST_DIR, "multiz45way-excerpt.maf") + +REPORTED_BLOCK_START = "20884758" + + +def run(args): + return subprocess.run([sys.executable, "-m", "mafutils"] + args, + check=False, cwd=REPO_ROOT, capture_output=True, text=True) + + +def maf_text(): + with open(MAF_FILE, "r", encoding="utf-8") as fp: + return fp.read() + + +def s_counts_per_block(text): + """Independent count: lines whose first field is 's', per 'a' block.""" + counts = [] + for line in text.splitlines(): + fields = line.split() + if not fields: + continue + if fields[0] == "a": + counts.append(0) + elif fields[0] == "s": + counts[-1] += 1 + return counts + + +def read_index(path): + with open(path, "r", encoding="utf-8") as fp: + lines = fp.read().splitlines() + return COMMON.readIndexHeader(path), [line.split("\t") for line in lines[1:]] + + +def build_index(maf_path): + result = run(["index", maf_path, "--progress-interval", "0"]) + assert result.returncode == 0, result.stdout + result.stderr + return COMMON.deriveBlockIndexPath(maf_path), COMMON.deriveScaffoldIndexPath(maf_path) + + +def copy_as(tmp_path, compression): + data = maf_text().encode("ascii") + if compression == "none": + path = os.path.join(str(tmp_path), "bird.maf") + with open(path, "wb") as fp: + fp.write(data) + elif compression == "gz": + path = os.path.join(str(tmp_path), "bird.maf.gz") + with gzip.open(path, "wb") as fp: + fp.write(data) + else: + path = os.path.join(str(tmp_path), "bird.maf.bgz") + with bgzf.BgzfWriter(path, "wb") as fp: + fp.write(data) + assert COMMON.detectCompression(path) == compression + return path + + +# --------------------------------------------------------------------------- +# Counts +# --------------------------------------------------------------------------- + + +def test_fixture_has_the_line_types_this_file_is_about(): + kinds = [line[:1] for line in maf_text().splitlines() if line] + assert kinds.count("s") == 72 and kinds.count("i") == 58 and kinds.count("e") == 12 + + +@pytest.mark.parametrize("compression", ["none", "gz", "bgzip"]) +def test_count_is_s_lines_only(tmp_path, compression): + block_idx, _ = build_index(copy_as(tmp_path, compression)) + header, rows = read_index(block_idx) + expected = s_counts_per_block(maf_text()) + + assert [int(r[5]) for r in rows] == expected + assert header["format"] == "3" + assert int(header["seq_lines"]) == sum(expected) == 72 + assert int(header["max_seqs"]) == max(expected) + assert int(header["blocks"]) == len(expected) == 14 + assert [int(r[5]) for r in rows if r[1] == REPORTED_BLOCK_START] == [2] + + +def test_offsets_bracket_blocks_with_that_many_s_lines(tmp_path): + maf = copy_as(tmp_path, "none") + block_idx, _ = build_index(maf) + _, rows = read_index(block_idx) + with open(maf, "rb") as fp: + data = fp.read() + for row in rows: + block = data[int(row[6]):int(row[7])] + assert block.startswith(b"a ") + assert sum(line.startswith(b"s") for line in block.split(b"\n")) == int(row[5]) + + +def test_comment_and_blank_lines_do_not_change_the_counts(tmp_path): + plain = copy_as(tmp_path, "none") + noisy = os.path.join(str(tmp_path), "noisy.maf") + out = [] + for line in maf_text().splitlines(keepends=True): + out.append(line) + if line.startswith("s "): + out.append("# a comment inside a block\n") + elif not line.strip(): + out.append("\n\n") + with open(noisy, "w", encoding="utf-8") as fp: + fp.writelines(out) + + _, plain_rows = read_index(build_index(plain)[0]) + _, noisy_rows = read_index(build_index(noisy)[0]) + # Offsets (columns 7-8) move with the inserted bytes; everything else must not. + assert [r[:6] for r in noisy_rows] == [r[:6] for r in plain_rows] + + +# --------------------------------------------------------------------------- +# Malformed MAFs: one clear error, and no index left behind +# --------------------------------------------------------------------------- + + +def _track_line_on_top(lines): + return lines[:1] + ["track name=birds\n"] + lines[1:] + + +def _indented_s_line(lines): + i = [n for n, line in enumerate(lines) if line.startswith("s ")][1] + return lines[:i] + [" " + lines[i]] + lines[i + 1:] + + +def _s_line_before_first_block(lines): + first_s = next(line for line in lines if line.startswith("s ")) + return lines[:1] + [first_s] + lines[1:] + + +def _block_without_s_lines(lines): + a_lines = [n for n, line in enumerate(lines) if line.startswith("a ")] + start, end = a_lines[2], a_lines[3] + return lines[:start] + [l for l in lines[start:end] if not l.startswith("s ")] + lines[end:] + + +@pytest.mark.parametrize("mutate, message", [ + (_track_line_on_top, "unrecognized line before the first block: 'track name=birds'"), + (_indented_s_line, "unrecognized line in block 1:"), + (_s_line_before_first_block, "'s' line before the first 'a' line"), + (_block_without_s_lines, "block 3 has no 's' line"), +], ids=["track-line", "indented-s-line", "s-before-first-a", "block-without-s"]) +def test_malformed_maf_errors_and_writes_no_index(tmp_path, mutate, message): + maf = os.path.join(str(tmp_path), "bad.maf") + with open(maf, "w", encoding="utf-8") as fp: + fp.writelines(mutate(maf_text().splitlines(keepends=True))) + + result = run(["index", maf, "--progress-interval", "0"]) + combined = result.stdout + result.stderr + assert result.returncode == 1, combined + assert message in combined + assert "No index was written" in combined + assert "Traceback" not in combined + assert not os.path.exists(COMMON.deriveBlockIndexPath(maf)) + assert not os.path.exists(COMMON.deriveScaffoldIndexPath(maf)) + + +# --------------------------------------------------------------------------- +# Older indexes (format 2 / headerless): refused only where the count is shown +# --------------------------------------------------------------------------- + + +def old_index_maf(tmp_path, kind): + """A MAF copy whose fresh index pair is rewritten as format 2, or headerless.""" + maf = copy_as(tmp_path, "none") + for path in build_index(maf): + with open(path, "r", encoding="utf-8") as fp: + lines = fp.readlines() + if kind == "format2": + assert " format=3 " in lines[0] + lines[0] = lines[0].replace(" format=3 ", " format=2 ") + else: + lines = lines[1:] + with open(path, "w", encoding="utf-8") as fp: + fp.writelines(lines) + return maf + + +@pytest.mark.parametrize("kind, described", [ + ("format2", "index format 2"), + ("headerless", "no index header"), +]) +def test_stats_refuses_an_old_index_unless_the_block_table_is_off(tmp_path, kind, described): + maf = old_index_maf(tmp_path, kind) + out = os.path.join(str(tmp_path), "out") + + refused = run(["stats", maf, "-o", out, "--progress-interval", "0"]) + combined = refused.stdout + refused.stderr + assert refused.returncode == 1, combined + assert described in combined and "--no-block-table" in combined + assert "Rebuild it with `mafutils index" in combined + assert not os.path.exists(out + ".overall.tsv") + + allowed = run(["stats", maf, "-o", out, "--no-block-table", "--progress-interval", "0"]) + assert allowed.returncode == 0, allowed.stdout + allowed.stderr + assert os.path.exists(out + ".overall.tsv") + + +@pytest.mark.parametrize("kind, described", [ + ("format2", "index format 2"), + ("headerless", "a headerless index"), +]) +def test_info_marks_old_counts_unavailable_and_exits_1(tmp_path, kind, described): + maf = old_index_maf(tmp_path, kind) + result = run(["info", maf, "--sample-blocks", "0"]) + assert result.returncode == 1, result.stdout + result.stderr + assert "blocks : 14" in result.stdout + assert f"seq lines : unavailable ({described} may count i/e/q lines)" in result.stdout + assert f"max seqs/block : unavailable ({described} may count i/e/q lines)" in result.stdout + assert "Rebuild it with `mafutils index" in result.stderr + + +def test_validate_reports_an_old_format_without_changing_the_verdict(tmp_path): + maf = old_index_maf(tmp_path, "format2") + result = run(["validate", maf]) + combined = result.stdout + result.stderr + assert result.returncode == 0, combined + assert "VERIFIED" in combined + assert "index format 2" in combined + + +def test_validate_reports_a_current_format(tmp_path): + maf = copy_as(tmp_path, "none") + build_index(maf) + result = run(["validate", maf]) + assert result.returncode == 0, result.stdout + result.stderr + assert "format : 3" in result.stdout + result.stderr + + +def test_fetch_and_gc_accept_a_format_2_index(tmp_path): + maf = old_index_maf(tmp_path, "format2") + bed = os.path.join(str(tmp_path), "region.bed") + with open(bed, "w", encoding="utf-8") as fp: + fp.write("chr6\t20884000\t20885000\tr1\n") + + fetched = run(["fetch", maf, bed, "-b", "id", "-o", os.path.join(str(tmp_path), "fetch"), + "--progress-interval", "0"]) + assert fetched.returncode == 0, fetched.stdout + fetched.stderr + assert os.path.isfile(os.path.join(str(tmp_path), "fetch", "r1.maf")) + + gc = run(["gc", maf, "-o", os.path.join(str(tmp_path), "gc"), "--progress-interval", "0"]) + assert gc.returncode == 0, gc.stdout + gc.stderr + + for result in (fetched, gc): + assert "i/e/q" not in result.stdout + result.stderr + + +# --------------------------------------------------------------------------- +# stats: with a format-3 index, its count and stats' own parse agree +# --------------------------------------------------------------------------- + + +def test_stats_index_and_parsed_counts_agree(tmp_path): + maf = copy_as(tmp_path, "none") + build_index(maf) + out = os.path.join(str(tmp_path), "out") + result = run(["stats", maf, "-o", out, "--progress-interval", "0"]) + assert result.returncode == 0, result.stdout + result.stderr + + with open(out + ".block.tsv", "r", encoding="utf-8") as fp: + rows = list(csv.DictReader(fp, delimiter="\t")) + assert len(rows) == 14 + assert [int(r["n_seq_lines_index"]) for r in rows] == [int(r["n_seq_lines_parsed"]) for r in rows] + assert sum(int(r["n_seq_lines_parsed"]) for r in rows) == 72 diff --git a/tests/test_info.py b/tests/test_info.py index 00544f9..e9abdd1 100644 --- a/tests/test_info.py +++ b/tests/test_info.py @@ -81,12 +81,13 @@ def test_validate_still_verifies_with_aggregates(tmp_path): 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.""" +def test_format_version_is_current(tmp_path): + """The aggregate keys are additive and did not bump the format (old readers + ignore unknown keys). Format 3 is the change in what the per-block + sequence count means: 's' lines only, no longer i/e/q lines too.""" maf = build_index(tmp_path) header = COMMON.readIndexHeader(COMMON.deriveBlockIndexPath(maf)) - assert header["format"] == "2" + assert header["format"] == "3" # --------------------------------------------------------------------------- diff --git a/tests/test_validate.py b/tests/test_validate.py index e8ca9a6..ee69681 100644 --- a/tests/test_validate.py +++ b/tests/test_validate.py @@ -52,7 +52,7 @@ def combined_output(result): def test_index_header_has_size_mtime_hash(): header = COMMON.readIndexHeader(INDEX_FILE) assert header is not None - assert header["format"] == "2" + assert header["format"] == "3" assert int(header["size"]) == os.path.getsize(MAF_FILE) assert header["hash"] == COMMON.computeFileHash(MAF_FILE) From 7c36a16fe64d039df8d4472812da52f340911969 Mon Sep 17 00:00:00 2001 From: gwct Date: Sun, 8 Nov 2026 21:08:51 -0500 Subject: [PATCH 3/3] Add possible future work to DEVELOPMENT.md - Duplicate checks in stats: assessed, but deferred over the temporary-disk cost on very large MAFs. - fetch copies i/e/q lines through trimmed blocks unchanged, so a q line no longer matches the length of its trimmed s line. - The e-line status could tell an indel from missing data. Co-Authored-By: Claude Opus 5.5 --- DEVELOPMENT.md | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index 74fd1b6..8b4988c 100644 --- a/DEVELOPMENT.md +++ b/DEVELOPMENT.md @@ -702,6 +702,55 @@ for user-facing installation/usage docs. backward-seek risk entirely rather than relying on any per-region processing order. +## Possible future work + +- **Duplicate checks in `stats`** (assessed 2026-11, not implemented). + "Duplicates" reported in `cactus-hal2maf --outType single` ("nodupes") MAFs + are expected, not a Cactus bug. That option guarantees at most one row per + genome per block, not that each non-reference base appears only once: where + the reference has several copies of a region, a single-copy segment of + another genome aligns at each of them. On a 27-genome hamster MAF this was + 2.2-3.7% of each non-reference genome's aligned bases, and 0 for the + reference. + - *A genome in more than one row of a block* is already counted per + species (`duplicated_blocks` in `species.tsv`, `n_duplicate_copies` in + `block.tsv`). Still missing, and cheap: + - an overall count of affected blocks (summing the species counts + double-counts); + - the first few locations; + - a log line worded for both nodupes and full MAFs. + - *Duplicate bases per genome* = aligned bases minus distinct bases, using + forward-strand coordinates (`srcSize - start - size` for `-` rows) and + keyed by (genome, sequence); optionally a BED of the duplicated intervals. + It needs every row's interval: + - merge intervals within each worker chunk (keeps ~36% of rows on a + hamster slice), counting overlaps within the chunk before merging; + - write compact binary records to one temporary file per genome; + - sort each genome in memory. + + Estimated at about 6 bytes of temporary disk per MAF row (~0.6 GB for the + 42.7 GB hamster MAF), so it should be opt-in. Check a large MAF's row + count first (`seq lines` in `mafutils info`). Deferred mainly over that + temporary-disk cost on very large MAFs. +- **`fetch` copies `i`/`e`/`q` lines through trimmed blocks unchanged** + (`trimMafBlock`, the `fields[0] != "s"` branch). In a block trimmed to a + region, a `q` line keeps its full-width quality string while its `s` line is + cut, so the two lengths no longer match. `i`-line base counts also go stale + at trimmed edges. `e` lines stay valid, and FASTA output never reads these + lines. Only MAFs that have them are affected: UCSC multiz MAFs, not cactus + ones. A fix would trim `q` lines with their `s` line and drop or recompute + `i` lines in trimmed blocks only, checked for byte-identity on whole blocks. +- **`e`-line status could tell an indel from missing data.** `e` marks a + species with no aligned bases in a block but bridged across it. Its status + says why: + - C: the flanking sequence is contiguous, so the region was deleted in that + species or inserted in the reference; + - I: unaligned bases in between; + - M: those bases are mostly N, i.e. missing data. + + `fetch` and `stats` treat all of these, and "not aligned at all", as the same + absence. This applies to UCSC MAFs only. + ## Development setup From a source checkout, an editable install picks up code changes without