diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index 192dfb1..8b4988c 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` @@ -155,6 +208,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`, @@ -183,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 @@ -296,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 @@ -599,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 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/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/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_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 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)