From b491cc363c0a1bbb62f1d66acbe3dec22b0734d3 Mon Sep 17 00:00:00 2001 From: gwct Date: Tue, 22 Sep 2026 10:20:01 -0400 Subject: [PATCH] Classify per-species loci in fetch and always emit the scaffold in FASTA headers --- DEVELOPMENT.md | 151 +++++++ README.md | 180 +++++++- mafutils/fetch.py | 517 +++++++++++++++++++++-- tests/loci-classes.maf | 97 +++++ tests/loci-classes.maf.block.idx | 7 + tests/loci-classes.maf.scaffold.idx | 2 + tests/test_loci.py | 613 ++++++++++++++++++++++++++++ 7 files changed, 1521 insertions(+), 46 deletions(-) create mode 100644 tests/loci-classes.maf create mode 100644 tests/loci-classes.maf.block.idx create mode 100644 tests/loci-classes.maf.scaffold.idx create mode 100644 tests/test_loci.py diff --git a/DEVELOPMENT.md b/DEVELOPMENT.md index 93e639b..2e921ef 100644 --- a/DEVELOPMENT.md +++ b/DEVELOPMENT.md @@ -296,6 +296,157 @@ for user-facing installation/usage docs. **warn and continue**. Real misuse is still caught loudly and earlier: index/MAF size+hash mismatch, a mismatched block/scaffold index pair, and per-scaffold extraction failures. +- **A stitched non-reference FASTA row is homology-by-alignment, not a called + locus — and its header coordinates are a bounding box.** This is the single + most confusing property of `fetch --fasta`, so: MAF blocks are contiguous in + the **reference** only. Each block is individually a proper alignment (verified + on 20,000 real blocks: every `s`-line's non-gap count equals its declared + size, every block rectangular), and a block boundary exists *precisely + because* contiguity broke somewhere — large indel, inversion, rearrangement, + assembly break, or a change in the aligned species set. Stitching a + multi-block region concatenates those independent alignments; **no padding is + inserted for the query's genomic distance**, so a 5 kb jump or a scaffold + change costs zero characters. Column correspondence is fully preserved (so it + is still an alignment, and phastCons/phyloP on it are sound), but each row is + not necessarily one genomic interval. + `writeFASTA` reports `min(start)`..`max(end)` over the contributing chunks, + which is therefore a bounding box. Observed on real data: a header claiming + `85525676-215275966` (129.75 Mb) for **89 emitted bases** — a 1.46-million-fold + overstatement. Measured over 2,000 real regions (25,599 region x non-reference + species pairs): 77.8% `single`, 10.7% `contiguous`, 10.6% `split` (median gap + 45 bp, max 335 Mb), 0.05% `multi_strand`, 0.8% `multi_scaffold`. So ~88.5% are + genuinely one locus and ~2.3% are not defensible as one. + `classifyChunks()` reports that classification and `--loci-table` emits the + per-chunk truth, so the indefensible cases can be filtered rather than + silently averaged into a span. **No gap threshold is hardcoded**, so the + number is reported and the judgment left to the analysis (same reasoning as + the advisory no-overlap thresholds). +- **A split gap is NOT an indel the alignment absorbed.** An earlier draft of + this note claimed a 45 bp split was "a small indel, still the same locus"; + the data disproves it. The aligner *does* thread query insertions inline as + reference-row gap columns, but only up to ~37 bp (156,928 such runs observed, + median 2, max 37, **zero at >=40**); past that it breaks the block instead. + Split gaps start at a median of 45 bp — immediately above that ceiling — so + they are sequence the aligner *declined to align*: real query bases silently + absent from the emitted FASTA while the bounding-box span still covers them. + Blocks tile the reference exactly (19,999/19,999 consecutive pairs adjacent, + zero gaps or overlaps), so a query jump at a boundary is purely query-side. + The consequence is exact, with zero exceptions across 27,599 records: + `single`/`contiguous` have header span == emitted bases and `split`/`multi_*` + never do — which is what makes the classification a usable filter. The + ~37 bp ceiling is an empirical property of this file's aligner and + parameters, not a MAF guarantee, so nothing hardcodes it. +- **Splits are mostly shared, not lineage-specific.** At a given block boundary + the number of species jumping is bimodal: 33.3% of boundaries have exactly + one species, but 17% have 13-14 of 14. Where several jump, ~30% agree on gap + size to within 10 bp (a real case: six species at 33/34/33/33/33/33 — one + ancestral indel inherited by descent). Nearly-everyone-jumps most + parsimoniously means a reference-specific deletion, which makes the + *element's* reference definition the questionable thing rather than one + species' copy. `summarizeElementLoci()` reports this as + `split.max.boundary.n.species` vs `n.split`. +- **Per-element classification lives in `maf_fetch_summary.tsv`, which is why + `trimMafBlock` returns chunks.** Classification is irreducibly per-species + (one real element was `contiguous` for 5 species, `split` for 6 and + `no_bases` for 3), but filtering a CNEE set needs a per-*element* verdict, and + that file is already one row per element. `trimMafBlock` was already computing + every trimmed `start`/`size`/`strand` to rewrite each `s` line and throwing + them away into the output text; returning them as a chunk list is what lets + the summary columns cover plain MAF output. Cost measured on 5,000 real + regions: **2.07s -> 2.21s (~6%)** in MAF mode, memory unchanged; FASTA mode + pays nothing new since it already accumulated chunks via `mafBlockToFasta`. + On real data 65.8% of elements are a clean single locus in all 15 species. + **One semantic difference to know:** for a MAF holding duplicate copies of a + species within one block, the MAF path counts every copy while the FASTA path + under `--fasta-dedupe most-seq` counts only the copy it kept. The production + file is `nodupes` (verified: zero repeated species across 20,000 blocks), so + the two agree there and the test suite asserts that equality. +- **`classifyChunks` excludes zero-size chunks, and `no_bases` exists because + `single` used to lie.** A species can carry an `s` line in a block whose + aligned sequence falls entirely outside the requested region once trimmed, + leaving an all-gap FASTA row and a zero-width header span. 5.46% of real loci + rows are zero-size, and **1,648 of them classified as a clean `single`** with + `max_gap=0` — so filtering `single|contiguous` to build a per-species BED + would have silently admitted ~1,650 zero-width intervals. Zero-size chunks + also invented phantom gaps: a real case (`Dicrostonyx_torquatus`, + `CM000994.3:10038193`) had a zero-size chunk in block 1 and 27 real bases in + block 3, which measured as a 718 bp `split` from a chunk containing no bases. + Such chunks are now dropped from the gap/class arithmetic but still emitted + as loci rows, so nothing is hidden. +- **Element-level and species-level filtering are very different yields; document + both.** The all-species filter (`n.split == n.multi.* == n.no.bases == 0`) + answers "is this element clean for everyone", which is what a concatenated + alignment needs -- 65.8% of elements on the 15-species rodent data. A + single-species product (e.g. a BED in one non-reference genome) should filter + that species' own rows in the loci table instead: 82-86% per species (mean + 83.7%) on the same data, and 98.6% reported on a more fragmented 45-way bird + alignment. Using the strict filter for a single-species product discards most + of the usable data, so the README shows both. +- **`interblock.distances` was removed; the check became a warning.** It + reported `next_ref_start - prev_ref_end` between consecutive blocks, which is + structurally zero on a normal MAF (19,999/19,999 consecutive pairs adjacent + on real data; the column read `NA`/`[0]`/`[0, 0]` and never anything else + across 2,000 real elements). Worse, it invited being read as a contiguity + check when it describes only the reference — the reference is contiguous by + construction while a species can jump kilobases at the same boundary, which + is the whole asymmetry the `split.*` columns exist to report. The computation + is retained and now raises a `ref-coverage-gap` warning when any distance is + non-zero, meaning the region holds reference positions no block covers (a + negative distance would mean blocks overlap in the reference, which should + not happen). Nothing is lost: `ref.block.bases` falls below `end - start` by + exactly the missing amount (verified: region width 9, `ref.block.bases` 3, + interblock distance `[6]`). +- **`n.overlapping.blocks`/`block.lengths` became `ref.n.overlapping.blocks`/ + `ref.block.bases`.** The prefix marks them as reference-side measurements so + they are not mistaken for per-species facts, and `block.lengths` was actively + misleading: the plural implied a list when it is a single sum. Note the + no-output accounting in `fetchByBatch` reads summary **field 4** + (`ref.n.overlapping.blocks`) positionally, which is unaffected by removing + field 6 — but it is positional, so check it if the leading columns ever move. +- Documentation examples that filter these files now look columns up **by + name** rather than by position. A positional example had already gone stale + once (`$13=="single"` for the loci table, after `block_index`/`gap_to_next` + shifted `class` to column 15), which is exactly the failure the name-based + form prevents. +- **Zero-size rows are a trap for consumers aggregating the loci table.** They + are emitted deliberately (nothing hidden) and excluded from `class`, so a + `single`-classified (element, species) can still carry a zero-size row -- + *possibly on a different scaffold*. A real case: `lod=251` / + `Dicrostonyx_torquatus` has a zero-size chunk on `JAQHUG010000114.1` and 106 + real bases on `JAQHUG010000001.1`, classified `single` (correctly). An awk + recipe that merged min/max across all rows without filtering `chunk_size > 0` + spanned both scaffolds and produced a 5.9 Mb interval for a ~250 bp element. + Any documented recipe must include the `chunk_size > 0` guard; the README's + does, and its output was verified against the reference element widths + (median 109 bp vs 129 bp reference, max 3,215 vs 3,213). +- **Coordinates are reported exactly as the MAF states them, never converted.** + For `strand = -`, MAF's `start` counts in the reverse-complemented source; the + old code had an `if strand == "+": ... else: ...` whose branches were + **identical** and whose comment claimed forward-chromosome reporting it never + did. The dead branch is gone. Forward conversion was considered and + deliberately rejected: the sequence must stay exactly what was aligned (that + is what makes the columns homologous, and it feeds molecular-evolution + analysis), and two coexisting coordinate frames would be worse than one. The + loci table emits `src_size` instead, so forward coordinates are one + subtraction away (`src_size - (start+size)`, `src_size - start`) without + mafutils picking a frame. +- **The FASTA dict key used to double as the header identity, which silently + dropped the scaffold.** `mafBlockToFasta` keyed records by `species` when + deduping (`use_species_keys = bool(expected_species) or fasta_dedupe != "none"`) + and by the full `species.scaffold` src otherwise — and `writeFASTA` built the + header from whatever the key was. So `-fh species-coords` emitted the scaffold + under `--fasta-dedupe none` but not under `most-seq`, making the same header + mode mean different things depending on unrelated flags, and making a + non-reference record unresolvable to a locus. `src`/`src_size` now travel + inside the record (swapped together with start/end/strand, so the copy dedupe + keeps is the one reported), the scaffold is built from the chunks via + `fastaScaffoldField()`, and multi-scaffold records comma-join their sources + rather than silently naming one. `species-only` stays bare on purpose — it + exists so headers match tree tip labels for phast/PhyloAcc. +- `writeFASTA` assigned `id_str` inside `if region["id"]` but used it + unconditionally, so `-fh species-coords-id` with a 3-column BED (no name + column) died with `UnboundLocalError`. The id suffix is now appended only when + present. - The vendored `bgzf.BgzfReader` has no `.name` attribute (unlike `gzip.GzipFile`), so `fetch.py` gets the MAF's display filename from the known `maf_file` path (`WORKER_MAF_FILE`), not from the open file handle. diff --git a/README.md b/README.md index 3d9b0e1..5f60237 100644 --- a/README.md +++ b/README.md @@ -177,17 +177,195 @@ Options: | `--basename`, `-b` | Output basename strategy: `id`, `coords`, or `count` | | `--output`, `-o` | Output directory or output filename in single-output mode | | `--fasta`, `-f` | Write FASTA instead of MAF | -| `--fasta-header`, `-fh` | FASTA header format: `species-coords-id`, `species-coords`, or `species-only` | +| `--fasta-header`, `-fh` | FASTA header format: `species-coords-id`, `species-coords`, or `species-only`. The coords-bearing modes always include the source scaffold (`species.scaffold`); `species-only` is deliberately bare so headers match tree tip labels | | `--expected-species` | Comma-separated expected species list for FASTA filling | | `--expected-species-file` | File with one expected species name per line | | `--fasta-dedupe` | FASTA duplicate handling: `none` or `most-seq` | | `--processes`, `-p` | Number of worker processes (see Compression above — plain gzip always runs single-process) | | `--mode`, `-m` | Fetch mode: `block` or `scaffold` | | `--scaffold-subdirs` | Group output files into subfolders named by reference scaffold (`//`) instead of one flat directory | +| `--loci-table` | Also write `/maf_fetch_loci.tsv`: one row per region/species/contributing-block with each species' source scaffold, coordinates, `srcSize`, and a locus classification (see *Per-species coordinates* below). Requires `--fasta` | | `--verbose` | Emit warning lines from each completed batch | | `--profile` | Log internal timing breakdowns | | `--verify-hash` | Verify the index's stored content hash against the MAF file (see Index Integrity above) | +### Reference-side columns in `maf_fetch_summary.tsv` + +Three of the original columns describe the **reference**, not any species, and +the `ref.` prefix marks them as such: + +| column | meaning | +|---|---| +| `ref.n.overlapping.blocks` | how many alignment blocks the region overlaps | +| `ref.block.bases` | reference bases in the region that alignment blocks cover | + +`ref.block.bases` normally equals `end - start`. It falls short exactly when +the region contains reference positions no block covers, and in that case +`fetch` emits a `ref-coverage-gap` warning naming the distances involved. + +This distinction matters: **the reference is contiguous by construction, so +these columns say nothing about whether a species' sequence is contiguous.** A +species can jump kilobases at a boundary where the reference does not. The old +`interblock.distances` column was removed for exactly this reason — it +reported reference-side gaps, which are all zeros on a normal MAF (measured: +19,999/19,999 consecutive block pairs adjacent), and reading it as a +contiguity check was misleading. Per-species contiguity lives in the `n.*` / +`split.*` columns below and in the loci table's `gap_to_next`. + +### Per-element locus classification (always written) + +`maf_fetch_summary.tsv` carries one row per BED element and reports how each +species' contribution to that element looks. These columns are always written +(no flag needed) and work with MAF or FASTA output: + +| column | meaning | +|---|---| +| `n.species` | species contributing to this element | +| `n.single`, `n.contiguous` | species whose contribution **is** one locus — for these, the header span equals the bases emitted | +| `n.split` | species interrupted by sequence the aligner declined to align | +| `n.multi.scaffold`, `n.multi.strand` | species stitched from different source sequences or strands | +| `n.no.bases` | species present in the alignment but contributing zero bases to this element (an all-gap row) | +| `split.bases.max` | worst species' total unaligned bases, i.e. how much sequence its row is missing | +| `split.gap.median` | median of every individual gap in the element | +| `split.boundaries` | distinct block boundaries with any split | +| `split.max.boundary.n.species` | most species splitting at a single boundary | +| `split.max.boundary.gap.spread` | `max - min` of those species' gaps at that same boundary | + +Both extent columns matter: one real element had a median gap of 63 bp across +13 split species but one species off by 335 Mb. The max alone would condemn an +element where 12 species are nearly fine; the median alone would wave through a +catastrophic misalignment. + +The concentration columns answer *"is one shared event responsible?"*. Splits +occur at block boundaries, which are single reference positions, so species +jumping at the same boundary are interrupted at the same place. Compare +`split.max.boundary.n.species` to `n.split`: equal, with a small +`gap.spread`, means one ancestral indel (a real case had six species at +33/34/33/33/33/33 — spread 1); much lower means independent lineage-specific +events. Where nearly every species jumps at one boundary, the most +parsimonious reading is a reference-specific deletion, which makes the +element's own reference definition the questionable thing. + +**Filter at the granularity your product actually needs** — these are very +different yields, and using the strict one for a single-species product +throws away most of your data. + +These examples look columns up by **name**, so they keep working if the +schema changes: + +*All species* — "safe to claim a locus in every species", e.g. for a +concatenated alignment you want uniformly clean: + +```bash +awk -F'\t' 'NR==1{for(i=1;i<=NF;i++)c[$i]=i; print; next} + $c["n.split"]==0 && $c["n.multi.scaffold"]==0 && + $c["n.multi.strand"]==0 && $c["n.no.bases"]==0' maf_fetch_summary.tsv +``` + +*One species* — for a per-species product such as a BED in a single +non-reference genome, filter that species' own rows in the loci table +instead; an element unusable for one species is usually fine for the rest. + +Two things to get right when aggregating those rows: + +1. **The table is one row per contributing block**, so a `contiguous` species + with two chunks yields two rows for one element. Since `contiguous` means + the chunks are adjacent by definition, merging them (min start, max end) + is lossless and gives one interval per element. On one real dataset a + species had 1,703 clean chunk rows across 1,462 elements — counting rows + would overstate the element count by 16%. +2. **Skip `chunk_size == 0` rows.** They are emitted so nothing is hidden, but + they contribute no sequence and are excluded from `class` — so a `single` + record can still have a zero-size row, *possibly on a different scaffold*. + Aggregating without filtering them produced a nonsense 5.9 Mb interval in + testing. + +A complete per-species BED, merged and converted to forward-strand +coordinates (verified: 1,462 intervals, median width 109 bp against a +reference element median of 129 bp): + +```bash +awk -F'\t' 'NR==1{for(i=1;i<=NF;i++)c[$i]=i; next} + $c["species"]=="YOUR_SPECIES" && $c["chunk_size"]+0>0 && + ($c["class"]=="single" || $c["class"]=="contiguous") { + k=$c["region_scaffold"]":"$c["region_start"] + s=$c["chunk_start"]+0; e=s+$c["chunk_size"] + if(!(k in lo) || shi[k]) hi[k]=e + scaf[k]=$c["src_scaffold"]; str[k]=$c["chunk_strand"] + sz[k]=$c["src_size"]+0; id[k]=$c["region_id"] + } + END{ for(k in lo){ + if(str[k]=="-"){ a=sz[k]-hi[k]; b=sz[k]-lo[k] } else { a=lo[k]; b=hi[k] } + print scaf[k]"\t"a"\t"b"\t"id[k]"\t0\t"str[k] + }}' maf_fetch_loci.tsv > species.bed +``` + +The `str[k]=="-"` branch is the forward-strand conversion described above — +`src_size - end` .. `src_size - start`, since mafutils reports MAF-frame +coordinates and leaves the choice of frame to you. + +On one real 15-species dataset the all-species filter kept **65.8%** of +elements while the per-species filter kept **82-86%** (mean 83.7%) depending +on the species; on a more fragmented 45-way alignment the per-species figure +was reported at **98.6%**. The gap is dataset-specific, so measure yours +rather than assuming either number. + +### Per-species coordinates (`--loci-table`) + +FASTA headers name the source scaffold, but for a **non-reference** species a +header's `start-end` is a *bounding box* over every block that contributed, not +necessarily one locus. MAF blocks are contiguous in the reference only: a +species' chunks can be separated, inverted, or on different scaffolds, and +stitching concatenates them without padding. On real data (2,000 regions, 15 +species) 88.5% of region/species pairs are unambiguously one locus, but ~11% +are split and ~0.8% span more than one scaffold. + +`--loci-table` writes the per-chunk truth so you can filter before treating a +row as a locus: + +```bash +mafutils fetch input.maf regions.bed -f -fh species-coords-id --loci-table -o out/ +``` + +`out/maf_fetch_loci.tsv` has one row per (region, species, contributing block): + +| column | meaning | +|---|---| +| `region_scaffold`, `region_start`, `region_end`, `region_id` | the requested reference region | +| `species`, `src_scaffold`, `src_size` | that species' source sequence for this chunk | +| `chunk_start`, `chunk_size`, `chunk_strand` | the chunk, exactly as the MAF states it | +| `chunk_index`, `n_chunks` | position among that species' chunks for this region | +| `block_index` | which block (ordinal within the element) this chunk came from | +| `gap_to_next` | signed gap to that species' next chunk; `.` for the last | +| `class` | `single`, `contiguous`, `split`, `multi_strand`, `multi_scaffold`, or `no_bases` | +| `max_gap` | largest signed gap between consecutive chunks; `0` for single/contiguous, `.` where undefined | + +`class` is constant within a (region, species) group, so filtering is a one-liner: + +```bash +awk -F'\t' 'NR==1{for(i=1;i<=NF;i++)c[$i]=i; next} + $c["class"]=="single" || $c["class"]=="contiguous"' out/maf_fetch_loci.tsv +``` + +**Coordinates are reported exactly as the MAF states them.** For `strand = -` +that means they count in the reverse-complemented source, which is MAF's own +convention. `src_size` is included so forward-genomic coordinates are one +subtraction away — `src_size - (chunk_start + chunk_size)` to +`src_size - chunk_start` — rather than mafutils choosing a frame for you. No +gap threshold is applied either. A `split` gap is **not** an indel the +alignment absorbed: the aligner threads query insertions inline as +reference-row gap columns only up to ~37 bp, and breaks the block past that, so +a gap between chunks is query sequence the aligner *declined to align* — real +bases absent from the FASTA while the header's span still covers them. The +consequence holds with no exceptions: `single`/`contiguous` records have +`header span == emitted bases`, `split`/`multi_*` never do. Real gaps range +from a median of 45 bp to 335 Mb, so `max_gap` is reported and the judgment is +yours. + +The table is roughly (regions x species) rows, so it is off by default — a +979k-region BED at 15 species produces ~13.7M rows. + On plain-gzip input, `fetch` precomputes the set of blocks needed across *all* regions, decodes each exactly once in strictly-ascending file order, then serves every region from that cache — this avoids the backward seeks diff --git a/mafutils/fetch.py b/mafutils/fetch.py index 8116148..5f7aaac 100644 --- a/mafutils/fetch.py +++ b/mafutils/fetch.py @@ -67,6 +67,8 @@ import logging import bisect import time +import statistics +import tempfile import traceback from enum import Enum from types import SimpleNamespace @@ -94,6 +96,7 @@ WORKER_OUTPUT = None WORKER_SINGLE_OUTPUT = None WORKER_SCAFFOLD_SUBDIRS = None +WORKER_LOCI_TMP_DIR = None WORKER_AS_FASTA = None WORKER_FASTA_HEADER = None WORKER_EXPECTED_SPECIES = None @@ -181,6 +184,7 @@ def initBatchWorker( scaffold_subdirs=False, index_file=None, scaffold_runs=None, + loci_tmp_dir=None, ): global WORKER_HEADER global WORKER_MAF_FILE @@ -191,6 +195,7 @@ def initBatchWorker( global WORKER_OUTPUT global WORKER_SINGLE_OUTPUT global WORKER_SCAFFOLD_SUBDIRS + global WORKER_LOCI_TMP_DIR global WORKER_AS_FASTA global WORKER_FASTA_HEADER global WORKER_EXPECTED_SPECIES @@ -212,6 +217,7 @@ def initBatchWorker( WORKER_OUTPUT = output WORKER_SINGLE_OUTPUT = single_output WORKER_SCAFFOLD_SUBDIRS = scaffold_subdirs + WORKER_LOCI_TMP_DIR = loci_tmp_dir WORKER_AS_FASTA = as_fasta WORKER_FASTA_HEADER = fasta_header WORKER_EXPECTED_SPECIES = expected_species @@ -497,12 +503,25 @@ def getMAFHeader(maf_file, maf_compression): def speciesFromSrc(src): return src.split(".", 1)[0] if "." in src else src + +def scaffoldFromSrc(src): + """ + The source-sequence part of a MAF "species.scaffold" src, or "" if the src + carries no scaffold at all. Mirrors gc.chromFromSrc, except that a src with + no "." yields "" rather than echoing the species back -- callers use this to + decide whether there is a scaffold to put in a header, and echoing the + species would produce "species.species". + """ + return src.split(".", 1)[1] if "." in src else "" + ############################################################################# WARNING_LABELS = { "missing-scaffold": "missing-scaffold", + "multi-scaffold": "multi-scaffold", "multiple-strands": "multiple-strands", "no-overlap": "no-overlap", + "ref-coverage-gap": "ref-coverage-gap", "trim-mismatch": "trim-mismatch", } @@ -728,24 +747,37 @@ def mafBlockToFasta(block_text, region, dedupe_mode="none", use_species_keys=Fal start = int(fields[2]) size = int(fields[3]) strand = fields[4] - #srcSize = fields[5] + src_size = int(fields[5]) seq = fields[6] block_len = len(seq) key = species if use_species_keys else src nongap = sum(1 for c in seq if c != "-") - - # Compute 1-based closed interval for clarity (as in FASTA tools) - if strand == "+": - end = start + size - else: - # For negative strand, coordinates can be reported as on "forward" chromosome (for clarity, as in MAF) - end = start + size + + # MAF-frame coordinates, reported exactly as the MAF states them: + # 0-based half-open [start, start+size). For strand "-" these count + # in the reverse-complemented source, which is MAF's convention -- + # deliberately NOT converted to forward-genomic here. src_size is + # carried alongside so a consumer that wants forward coordinates can + # compute them (src_size - (start+size), src_size - start) without + # mafutils baking in that interpretation. + end = start + size + + # scaffold/src_size travel WITH start/end/strand so that whichever + # copy dedupe keeps, its own source is the one reported. Previously + # the scaffold lived only in the dict key, so keying by species + # (--fasta-dedupe most-seq, or any --expected-species) discarded it + # entirely and the header could not name a locus. + record = { + 'seq': seq, 'src': src, 'scaffold': scaffoldFromSrc(src), + 'start': start, 'end': end, 'size': size, 'strand': strand, + 'src_size': src_size, 'nongap': nongap, + } if key in fasta_lines and dedupe_mode == "most-seq": if nongap > fasta_lines[key]["nongap"]: - fasta_lines[key] = {'seq': seq, 'start': start, 'end': end, 'strand': strand, 'nongap': nongap} + fasta_lines[key] = record elif key not in fasta_lines: - fasta_lines[key] = {'seq': seq, 'start': start, 'end': end, 'strand': strand, 'nongap': nongap} + fasta_lines[key] = record elif dedupe_mode == "none": # Preserve previous behavior when dedupe is off by keeping src keys unique. # If duplicate src appears, keep first occurrence. @@ -759,40 +791,302 @@ def mafBlockToFasta(block_text, region, dedupe_mode="none", use_species_keys=Fal ############################################################################# -def writeFASTA(fasta_seqs, region, fasta_stream, BATCHLOG, fasta_header=False, verbose=False, warning_state=None): +def 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 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) + ] + + +ELEMENT_LOCI_HEADERS = [ + "n.species", "n.single", "n.contiguous", "n.split", + "n.multi.scaffold", "n.multi.strand", "n.no.bases", + "split.bases.max", "split.gap.median", + "split.boundaries", "split.max.boundary.n.species", + "split.max.boundary.gap.spread", +] + + +def 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 = classifyChunks(chunks) + if locus_class in counts: + counts[locus_class] += 1 + + gaps = 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, + ] + + +def fastaScaffoldField(chunks): + """ + The scaffold portion of a FASTA header, built from the chunks rather than + from the dict key. Distinct scaffolds are comma-joined in block order, so a + record stitched from two source sequences says so instead of silently + naming one of them. Returns "" when no chunk carries a scaffold. + """ + seen = [] + for c in chunks: + scaf = c.get('scaffold') or "" + if scaf and scaf not in seen: + seen.append(scaf) + return ",".join(seen) + + +def writeFASTA(fasta_seqs, region, fasta_stream, BATCHLOG, fasta_header=False, verbose=False, warning_state=None, loci_stream=None): - fasta_output = {} for sp, details in fasta_seqs.items(): + chunks = details.get('chunks') or [] + species = speciesFromSrc(sp) + locus_class, max_gap = classifyChunks(chunks) + if fasta_header == "species-only": - header = f">{speciesFromSrc(sp)}" + # Deliberately bare: this mode exists so headers match tree tip + # labels for phast/PhyloAcc, so it never gains a scaffold. + header = f">{species}" else: - if details['starts'] and details['ends']: - region_start = min(details['starts']) - region_end = max(details['ends']) + if chunks: + region_start = min(c['start'] for c in chunks) + region_end = max(c['end'] for c in chunks) else: region_start = region['start'] region_end = region['end'] - strand = list(set(details['strands'])) if details['strands'] else ["."] - - if len(strand) == 1: - region_strand = strand[0] + + strands = {c['strand'] for c in chunks if c['strand'] is not None} + if len(strands) == 1: + region_strand = next(iter(strands)) else: region_strand = "." + if len(strands) > 1: + appendWarning( + warning_state, + "multiple-strands", + f"Multiple strands found for species {sp} in region {region['scaffold']}:{region['start']}-{region['end']}. Using '.' in header.", + ) + + # The scaffold comes from the chunks, not the dict key: keying by + # species (--fasta-dedupe most-seq / --expected-species) used to + # drop it entirely, so the header could not name a locus at all. + scaffold_field = fastaScaffoldField(chunks) + name = f"{species}.{scaffold_field}" if scaffold_field else species + + if locus_class == "multi_scaffold": appendWarning( warning_state, - "multiple-strands", - f"Multiple strands found for species {sp} in region {region['scaffold']}:{region['start']}-{region['end']}. Using '.' in header.", + "multi-scaffold", + f"Species {species} in region {region['scaffold']}:{region['start']}-{region['end']} " + f"is stitched from {len(set(c['scaffold'] for c in chunks))} source sequences " + f"({scaffold_field}); its header coordinates are a bounding box, not a locus.", ) + # start-end is min..max over the contributing chunks -- a bounding + # box. It equals the emitted sequence's extent only when the chunks + # are contiguous (see classifyChunks); the loci table carries the + # per-chunk truth. if fasta_header == "species-coords": - header = f">{sp}:{region_start}-{region_end}({region_strand})" - elif fasta_header == "species-coords-id": - if "id" in region and region["id"]: - id_str = f"id:{region['id']}" - header = f">{sp}:{region_start}-{region_end}({region_strand}) {id_str}" + header = f">{name}:{region_start}-{region_end}({region_strand})" + else: # species-coords-id + header = f">{name}:{region_start}-{region_end}({region_strand})" + # A BED with no name column leaves region["id"] as None; this + # used to leave id_str unbound and crash with UnboundLocalError. + if region.get("id"): + header += f" id:{region['id']}" fasta_stream.write(f"{header}\n{''.join(details['seq'])}\n") + if loci_stream is not None: + writeLociRows(loci_stream, region, species, chunks, locus_class, max_gap) + + +LOCI_TABLE_HEADERS = [ + "region_scaffold", "region_start", "region_end", "region_id", + "species", "src_scaffold", "src_size", + "chunk_start", "chunk_size", "chunk_strand", + "chunk_index", "n_chunks", "block_index", "gap_to_next", + "class", "max_gap", +] + + +def writeLociRows(loci_stream, region, species, chunks, locus_class, max_gap): + """ + One row per (region, species, chunk), with the per-(region, species) class + denormalized onto each row so a single file serves both filtering and + per-species BED construction. + + Coordinates are MAF-frame, exactly as the MAF states them; src_size is + included so forward-genomic coordinates are one subtraction away + (src_size - (start+size), src_size - start) without mafutils choosing a + frame on the caller's behalf. + + block_index and gap_to_next make boundary identity explicit: a split + happens AT a block boundary, which is one reference position, so two + species sharing a block_index are interrupted at the same place. Together + they make per-boundary species counts and gap agreement derivable + downstream instead of frozen into the tool. Rows are ordered by + chunk_start, so a non-ascending block_index itself flags a chunk running + against the reference (7 of 3,744 real gaps). + """ + gap_str = "." if max_gap is None else str(max_gap) + region_id = region.get("id") or "." + ordered = sorted(chunks, key=lambda c: c['start']) + gaps = chunkGaps(chunks) + n = len(ordered) + # chunkGaps only covers contributing (size>0) chunks on one scaffold/strand, + # so it is shorter than `ordered` whenever zero-size or heterogeneous chunks + # are present; those rows get "." rather than a misleading number. + contributing = [c for c in ordered if c['size'] > 0] + gap_by_start = {} + if len(gaps) == len(contributing) - 1: + for i, gap in enumerate(gaps): + gap_by_start[contributing[i]['start']] = gap + + for i, c in enumerate(ordered, start=1): + gap_to_next = gap_by_start.get(c['start']) + loci_stream.write("\t".join([ + region["scaffold"], str(region["start"]), str(region["end"]), region_id, + species, c.get('scaffold') or ".", + "." if c.get('src_size') is None else str(c['src_size']), + str(c['start']), str(c['size']), c.get('strand') or ".", + str(i), str(n), + str(c.get('block_index', i)), + "." if gap_to_next is None else str(gap_to_next), + locus_class, gap_str, + ]) + "\n") + + ############################################################################# def trimMafBlock(block_text, bed_start, bed_end, BATCHLOG, verbose=False, warning_state=None): @@ -804,9 +1098,15 @@ def trimMafBlock(block_text, bed_start, bed_end, BATCHLOG, verbose=False, warnin The MAF "s" line format: s - Returns the trimmed block as a string, or None if no overlap occurs. + Returns (trimmed_block_text, ref_info, chunks), or None if no overlap + occurs. `chunks` is one record per trimmed "s" line (species, scaffold, + trimmed start/size, strand, srcSize) -- every value already computed here + to rewrite the line, so collecting them costs no extra parsing. That is + what lets the per-element locus summary cover plain MAF output as well as + FASTA, which would otherwise need a second parse of every block. """ + chunks = [] lines = block_text.splitlines() if len(lines) < 2: return None # Invalid block @@ -886,6 +1186,15 @@ def trimMafBlock(block_text, bed_start, bed_end, BATCHLOG, verbose=False, warnin new_fields[3] = str(new_ref_size) new_fields[6] = ref_seq[ref_col_start:ref_col_end] trimmed_lines.append(" ".join(new_fields)) + chunks.append({ + 'species': speciesFromSrc(fields[1]), + 'scaffold': scaffoldFromSrc(fields[1]), + 'start': new_ref_start, + 'size': new_ref_size, + 'end': new_ref_start + new_ref_size, + 'strand': fields[4], + 'src_size': int(fields[5]) if len(fields) > 5 and fields[5].isdigit() else None, + }) else: species_seq = fields[6] @@ -908,10 +1217,19 @@ def trimMafBlock(block_text, bed_start, bed_end, BATCHLOG, verbose=False, warnin new_fields[3] = str(trimmed_size) new_fields[6] = trimmed_seq trimmed_lines.append(' '.join(new_fields)) + chunks.append({ + 'species': speciesFromSrc(fields[1]), + 'scaffold': scaffoldFromSrc(fields[1]), + 'start': new_start, + 'size': trimmed_size, + 'end': new_start + trimmed_size, + 'strand': strand, + 'src_size': int(fields[5]) if len(fields) > 5 and fields[5].isdigit() else None, + }) #print(f"DEBUG: bed_interval=[{bed_start},{bed_end}), block_ref=[{block_ref_start},{block_ref_start + block_ref_size}), extracted_cols=[{ref_col_start}:{ref_col_end}], ref_seq_sample='{ref_seq[ref_col_start:ref_col_end] if ref_col_start is not None else 'None'}'", file=sys.stderr) - return "\n".join(trimmed_lines), [new_ref_start, new_ref_end, new_ref_size] + return "\n".join(trimmed_lines), [new_ref_start, new_ref_end, new_ref_size], chunks ############################################################################# @@ -945,6 +1263,7 @@ def fetchByRegion( warning_state=None, profile_state=None, scaffold_subdirs=False, + loci_stream=None, ): """ Worker function to process a single BED region: @@ -966,6 +1285,9 @@ 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) + block_ordinal = 0 total_ref_bases = 0 # Total ref bases across all blocks if as_fasta: @@ -993,7 +1315,9 @@ def fetchByRegion( "missing-scaffold", f"{scaffold}:{bed_start}-{bed_end} No index entries for scaffold.", ) - summary = f"{region_str}\t0\t0\tNA\t0" + # Same column count as the normal path: no blocks means no species and + # no splits, so every element-locus column is zero. + summary = f"{region_str}\t0\t0\t0\t" + "\t".join(["0"] * len(ELEMENT_LOCI_HEADERS)) return summary if not single_output else [] blocks_written = 0 @@ -1030,7 +1354,23 @@ def fetchByRegion( addProfile(profile_state, "time_trim", time.perf_counter() - trim_timer_start) if trimmed_result is None: continue - trimmed, trimmed_info = trimmed_result + trimmed, trimmed_info, block_chunks = trimmed_result + + # Accumulate per-species chunks for EVERY output format. On the FASTA + # path the authoritative records come from mafBlockToFasta below, which + # is dedupe-aware (the chunk must describe the copy --fasta-dedupe + # most-seq actually kept); here they come straight from the trimmer, so + # MAF output gets the same per-element locus summary at no extra + # parsing cost. NOTE the one semantic difference: for a MAF holding + # duplicate copies of a species in one block, this counts every copy + # 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. + for chunk in block_chunks: + chunk['block_index'] = block_ordinal + element_chunks[chunk.pop('species')].append(chunk) if not single_output and trimmed is not None: block_stats.append(trimmed_info) @@ -1076,9 +1416,11 @@ def fetchByRegion( if sp not in fasta_seqs: # Backfill for all previous blocks fasta_seqs[sp]['seq'] = [getFillString(fill_cache, "-", l) for l in block_lengths[:-1]] - fasta_seqs[sp]['starts'] = [] - fasta_seqs[sp]['ends'] = [] - fasta_seqs[sp]['strands'] = [] + # One record per contributing block, in block order, as + # the single source of truth: the header's coordinates, + # its scaffold field, and the locus classification are + # all derived from this rather than from parallel lists. + fasta_seqs[sp]['chunks'] = [] # After establishing all species, append current block or pad as needed for sp in species_order: @@ -1088,13 +1430,22 @@ def fetchByRegion( fill_char = "N" if sp in expected_species_set else "-" seq = getFillString(fill_cache, fill_char, block_len) fasta_seqs[sp]['seq'].append(seq) - if sp in return_fasta_seqs: - if return_fasta_seqs[sp]['start'] is not None: - fasta_seqs[sp]['starts'].append(return_fasta_seqs[sp]['start']) - if return_fasta_seqs[sp]['end'] is not None: - fasta_seqs[sp]['ends'].append(return_fasta_seqs[sp]['end']) - if return_fasta_seqs[sp]['strand'] is not None: - fasta_seqs[sp]['strands'].append(return_fasta_seqs[sp]['strand']) + rec = return_fasta_seqs.get(sp) + # start is None for an --expected-species row that this + # block filled with N's: it contributed no real sequence, so + # it contributes no chunk. + if rec is not None and rec.get('start') is not None: + chunk = { + 'scaffold': rec.get('scaffold', ''), + 'start': rec['start'], + 'end': rec['end'], + 'size': rec.get('size', rec['end'] - rec['start']), + 'strand': rec.get('strand'), + 'src_size': rec.get('src_size'), + 'block_index': block_ordinal, + } + fasta_seqs[sp]['chunks'].append(chunk) + element_chunks[speciesFromSrc(sp)].append(chunk) if profile_state is not None: addProfile(profile_state, "time_fasta_stitch", time.perf_counter() - stitch_timer_start) else: @@ -1124,6 +1475,7 @@ def fetchByRegion( fasta_header=fasta_header, verbose=verbose, warning_state=warning_state, + loci_stream=loci_stream, ) if profile_state is not None: addProfile(profile_state, "time_write_fasta", time.perf_counter() - write_timer_start) @@ -1142,19 +1494,37 @@ def fetchByRegion( # Collect region info block_stats.sort() + # Consecutive blocks are expected to tile the reference exactly + # (verified on real data: 19,999/19,999 consecutive pairs adjacent), so + # this was reported as an all-zeros "interblock.distances" column that + # carried no information. It is a warning instead: a non-zero distance + # means the region contains reference bases with no alignment block, so + # ref.block.bases falls short of the region width and the extracted + # sequence is missing reference positions too. A negative distance would + # mean blocks overlap in the reference, which should not happen at all. interblock_spaces = [ block_stats[i][0] - block_stats[i-1][1] for i in range(1, len(block_stats)) ] if len(block_stats) > 1 else [] + nonzero_spaces = [d for d in interblock_spaces if d != 0] + if nonzero_spaces: + appendWarning( + warning_state, + "ref-coverage-gap", + f"{scaffold}:{bed_start}-{bed_end} blocks do not tile the reference " + f"(inter-block distances {nonzero_spaces}); the region contains reference " + f"bases with no alignment block.", + ) + num_blocks = len(block_stats) blocks_length = sum(info[2] for info in block_stats) - space_str = str(interblock_spaces) if interblock_spaces else "NA" if profile_state is not None: 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) - summary = (f"{region_str}\t{num_blocks}\t{blocks_length}\t{space_str}\t{len(species_order) if as_fasta else 'NA'}") + element_cols = "\t".join(str(v) for v in summarizeElementLoci(element_chunks)) + summary = (f"{region_str}\t{num_blocks}\t{blocks_length}\t{len(species_order) if as_fasta else 'NA'}\t{element_cols}") return summary else: @@ -1195,6 +1565,15 @@ def fetchByBatch( # to, by scaffold *name*) would break that: index scaffold order is MAF file # order, not alphabetical. index_fp = open(WORKER_INDEX_FILE, "r", encoding="utf-8") + # One loci file per BATCH, written straight to disk. Returning ~14 rows per + # region to the parent (as the region summary does) would not scale: a + # 979k-region run at 15 species is ~13.7M rows. The parent concatenates + # these in batch order -- see run_fetch. + loci_path = None + loci_fp = None + if WORKER_LOCI_TMP_DIR is not None: + loci_path = os.path.join(WORKER_LOCI_TMP_DIR, f"loci.task{batch_num}.tsv") + loci_fp = open(loci_path, "w", encoding="utf-8") scaffold_to_runs, runs, run_offsets = WORKER_SCAFFOLD_RUNS # Per-run resume points, so consecutive regions on a scaffold continue # forward through the index instead of rescanning from the run's start. @@ -1241,6 +1620,7 @@ def fetchByBatch( warning_state=warning_state, profile_state=profile_state, scaffold_subdirs=WORKER_SCAFFOLD_SUBDIRS, + loci_stream=loci_fp, ) if scan_state and "resume" in scan_state: @@ -1252,7 +1632,7 @@ def fetchByBatch( batch_results[result_idx] = region_result if not WORKER_SINGLE_OUTPUT: - # Field 4 is n.overlapping.blocks (see summary_headers). This + # Field 4 is ref.n.overlapping.blocks (see summary_headers). This # used to read field 1, which is the region's START coordinate # -- so every region beginning at coordinate 0 was counted as # "no output" despite having been extracted fine, inflating @@ -1269,6 +1649,8 @@ def fetchByBatch( pass finally: index_fp.close() + if loci_fp is not None: + loci_fp.close() if profile_state is not None: profile_state["time_batch_total"] = time.perf_counter() - batch_timer_start @@ -1279,6 +1661,7 @@ def fetchByBatch( "zero_overlap_regions": no_output_regions, "no_output_regions": no_output_regions, "batch_num": batch_num, + "loci_path": loci_path, "written_regions": written_regions, "warning_count": warning_state["count"], "warning_counts": warning_state["counts"], @@ -1372,6 +1755,13 @@ def run_fetch(args, cmdline="mafutils fetch"): sys.exit(1) # Validate mode + if args.loci_table and not args.fasta: + # The per-chunk records it reports are accumulated only on the FASTA + # stitching path; MAF output already carries species.scaffold in its + # s-lines, so nothing is missing there to reconstruct. + LOG.error("--loci-table requires --fasta (MAF output already carries species.scaffold in its s-lines).") + sys.exit(1) + if args.scaffold_subdirs and args.single_output: LOG.error("--scaffold-subdirs cannot be used with --single-output.") sys.exit(1) @@ -1677,6 +2067,9 @@ def run_fetch(args, cmdline="mafutils fetch"): missing_scaffold_regions = 0 profile_totals = initProfileState() if args.profile else None profiled_batches = 0 + loci_tmp = tempfile.TemporaryDirectory(prefix="maf_fetch_loci_", dir=args.output) if args.loci_table else None + loci_tmp_dir = loci_tmp.name if loci_tmp is not None else None + loci_paths = {} with ProcessPoolExecutor( max_workers=effective_processes, initializer=initBatchWorker, @@ -1697,9 +2090,17 @@ def run_fetch(args, cmdline="mafutils fetch"): args.scaffold_subdirs, index_file, scaffold_runs, + loci_tmp_dir, ), ) as executor, open(info_outfile, "w") as info_out: - summary_headers = ["scaffold", "start", "end", "basename", "n.overlapping.blocks", "block.lengths", "interblock.distances", "n.sequences"] + # ref.* prefix marks the columns measured on the REFERENCE, so they are + # not mistaken for per-species facts: the reference is contiguous by + # construction while a species can jump kilobases at the same boundary + # (that asymmetry is what the n.*/split.* columns exist to report). + # ref.block.bases was "block.lengths", whose plural implied a list when + # it is a single sum; it falls below (end - start) exactly when the + # region has reference bases no block covers. + summary_headers = ["scaffold", "start", "end", "basename", "ref.n.overlapping.blocks", "ref.block.bases", "n.sequences"] + ELEMENT_LOCI_HEADERS info_out.write("\t".join(summary_headers) + "\n") futures = {} @@ -1721,6 +2122,9 @@ def run_fetch(args, cmdline="mafutils fetch"): for key, value in result["profile"].items(): profile_totals[key] += value + if result.get("loci_path"): + loci_paths[result["batch_num"]] = result["loci_path"] + if not args.single_output: for r in result["results"]: info_out.write(r + "\n") @@ -1802,6 +2206,27 @@ def run_fetch(args, cmdline="mafutils fetch"): profile_totals["time_write_fasta"], ) + if loci_tmp is not None: + loci_outfile = os.path.join(args.output, "maf_fetch_loci.tsv") + rows = 0 + with open(loci_outfile, "w", encoding="utf-8") as loci_out: + loci_out.write("\t".join(LOCI_TABLE_HEADERS) + "\n") + # Concatenate in BATCH-NUMBER order, not filename order: sorting the + # paths as strings would put task10 before task2 (the ordering bug + # just fixed in stats' writeBlockTable). Batches partition the + # index-ordered regions, so this keeps the table deterministic and + # in genomic order. + for batch_num in sorted(loci_paths): + path = loci_paths[batch_num] + if not os.path.isfile(path): + continue + with open(path, "r", encoding="utf-8") as fp: + for line in fp: + loci_out.write(line) + rows += 1 + loci_tmp.cleanup() + LOG.info(f"Wrote {rows} loci row(s) to {loci_outfile}") + if not args.single_output and missing_scaffold_regions > 0: # Advisory, like the other non-overlap outcomes: fetch is normally # pointed at many intervals and a BED may legitimately carry scaffolds @@ -1884,6 +2309,7 @@ def fetch_command( processes: Annotated[int, typer.Option("--processes", "-p", help="Number of parallel processes to use (default: 1)")] = 1, mode: Annotated[FetchMode, typer.Option("--mode", "-m", help="Mode: 'block' to trim by regions, 'scaffold' to extract whole scaffolds (default: block)")] = FetchMode.block, scaffold_subdirs: Annotated[bool, typer.Option("--scaffold-subdirs", help="Group output files into subfolders named by reference scaffold (//) instead of one flat directory.")] = False, + loci_table: Annotated[bool, typer.Option("--loci-table", help="Also write /maf_fetch_loci.tsv: one row per region/species/contributing-block, with each species' source scaffold, MAF-frame coordinates, srcSize, and a locus classification (single/contiguous/split/multi_strand/multi_scaffold). Requires --fasta. Large: roughly (regions x species) rows.")] = False, single_output: Annotated[bool, typer.Option("--single-output", hidden=True)] = False, verbose: Annotated[bool, typer.Option("--verbose", help="Print every warning line after each batch completes.")] = False, profile: Annotated[bool, typer.Option("--profile", help="Log aggregate timing breakdowns for internal fetch steps.")] = False, @@ -1903,6 +2329,7 @@ def fetch_command( processes=processes, mode=mode.value, scaffold_subdirs=scaffold_subdirs, + loci_table=loci_table, single_output=single_output, verbose=verbose, profile=profile, diff --git a/tests/loci-classes.maf b/tests/loci-classes.maf new file mode 100644 index 0000000..4da1c99 --- /dev/null +++ b/tests/loci-classes.maf @@ -0,0 +1,97 @@ +##maf version=1 +a +s Mus_musculus.CM000994.3 3541841 56 + 195154279 TTGCTCCTTTTGTAAAACAGAAT---TTTCATTGCAAAACAAAAGAATAGATAATTCTT +s Rattus_norvegicus.CM070395.1 151308217 59 - 172190305 TTGCCCTCTTTGTAAAGCAGAATGATTTTGATTGCAAAACAGAACGATAGGTAATGAGT +s Phodopus_campbelli.chr5 215194483 4 - 232406260 TTG-----C-------------------------------------------------- +s Phodopus_sungorus_new.chr5 215606660 4 + 232644518 TTG-----T-------------------------------------------------- +s Phodopus_roborovskii.JAEOAF010000015.1 16898655 4 - 46562299 TTG-----C-------------------------------------------------- +s Cricetus_cricetus.OZ199625.1 165080875 4 + 494177253 ttg-----t-------------------------------------------------- +s Cricetulus_griseus.NC_048595.1 156951061 4 + 461620117 ttg-----t-------------------------------------------------- +s Dicrostonyx_torquatus.JAQHUG010000244.1 1627761 3 + 3161476 ttg-------------------------------------------------------- +s Arvicola_amphibius.NC_052057.2 11444849 3 - 123987292 TTG-------------------------------------------------------- +s Microtus_pennsylvanicus.CM073169.1 125811142 3 + 145900345 TTG-------------------------------------------------------- +s Neotoma_floridana.CM113784.1 121200281 47 + 180930968 ---------ATAGGAAGTGGAGT---TCTGTTCATAGGGATGCACCATGTGGAATGATT +s Onychomys_torridus.LR877189.1 121984033 47 - 167895514 ---------GTGGGAAGTGGAGT---TCTATTTGTAGGGATGCACCATATGGAGTGATT +s Peromyscus_californicus.VALE03000005.1 122460984 47 + 169319897 ---------GTGGGAAGTGAAGT---TATGTTCATAGGGATGCACCATATGGAGTGATT +s Peromyscus_maniculatus_North.CM113196.1 38376802 45 + 173669867 ----------TGGGGCGTGTAGT----CTGTTAGTAGGGAAGCACCATGTGGAATGATT + +a +s Mus_musculus.CM000994.3 3541897 163 + 195154279 TGTT--TGTCTTGAAATTTC-AAAAACAGCAATTATTGTGACT-AATATTTACTAAAATAATGTAAATCAAAC---ATCTTACTACCTTTAAAACTAAGAAA---AA--TG-ACTCATTTCT-AT-AGTCAGACTCTCAA-AATCTGTATGGTTCAAAGGATTGAGTAGTAAGTAGAAA +s Rattus_norvegicus.CM070395.1 151308276 167 - 172190305 GGTTTCTATCCTGAATTTTCAAAAAACAGCAATTATTGTGTCT-AACATTTACCAAAATAATGTAAATCAAACTT-ATCTTGTTACTTTAAAGG-TAAGAAA---AA--TG-GCTCATTTCT-AT-AGTCAGACTCTCAA-TATCTATAGGGTTCAAAGAATTGAGTAATAAGTAGAAA +s Phodopus_campbelli.chr5 215194525 168 - 232406260 TGTTTCTGTCTTGGAATTTC-AAAAACAGCAATTAATGTGTTT-AATATTTACCAAAATAGTGTATGTCAAATTT-ATCTTATTTCTTTAAAAACTAAGAAAGAAAA--AAAGCTCATTTCT-AT-AGCAGGACTCTCAG-AACCTGTAGGGTTCAAAGGACTGCACAGAAAATGG--- +s Phodopus_sungorus_new.chr5 215606702 167 + 232644518 TGTTTCTGCCTTGGAATTTC-AAAAACAGCAATTAATGTGTTT-AATATTTACCAAAATAGTGTATGTCAAATTT-ATCTTATTCCTTTAAAAACTAAGAAAGA-AA--AAAGCTCATTTCT-AT-AGCAGGACTCTCAG-AACCTGTAGGGTTCAAAGGACTGCATAGAAAATGG--- +s Phodopus_roborovskii.JAEOAF010000015.1 16898697 165 - 46562299 TGTTTCTGTCTTGGAATTTC-AAAAACAGCAATTAATGTGTTT-AACATTTATCAAAATAGAGTATGTCAACTTT-ATCTTATTCCTTTAAAAACCAAGAAA---AA--AAAGCTCATTTCT-AT-AGCAAGACTCTCAG-AAACTTTAGGGTTCAAAGGACTGAGTAGCAAATGG--- +s Cricetus_cricetus.OZ199625.1 124237172 164 + 494177253 TGGTGCTGTCTTGAAATTTC-AAAAATAGCAATTGTTGTGTCTTAATAATTACCAAAATAGTATAAGTTGAATTT-ATCTTACTACTTTAAAAAATAAGAAA---AT--TT--CCCAACTTC-TA-AATGGGACTCTCCA-AATATCTAGGGTTCAAAGGAAGGACTAGTAAGTAG--- +s Cricetulus_griseus.NC_048595.1 118041490 164 + 461620117 TGGTGCTGTCTTGAAATTTG-AAAAATAGCAATTGATGTGTCTTAATAATTACCAAAATAGTATAAGTTAAATTT-ATCTTACTACTTTAAAAACTAAGAAA---AA--TG--CCCAACTTC-TA-AATGGGAATCTCAA-AATCTCTAGGGTTCAAAGGAAGGACTAGTAAGTAG--- +s Dicrostonyx_torquatus.JAQHUG010000244.1 1628745 163 + 3161476 TGTTTCtgtcttgaaaattc--aaaacagcaatgaatgtgttt-aatatttatcaagatagtgtaactcaaatTTTatcttactattttaaaaactaaaaaa---at--GG--CTCATTTCT-AT-AATAGGACTCTCCA-AATCTGTAGGCTTCAAAGAACTGAATAATAAATAG--- +s Arvicola_amphibius.NC_052057.2 11445666 163 - 123987292 tattTCtgtcttgaaaattc-aaaaacaacaattaatgtgttt-aatatttatcaaactagtgtaaggcaaatTTTatcttactattttaaaaactaaaaag---aa--gg--ctCATTTCT-AT-AATTGGACTCTCC--AATCTGTTGGTTTTAAAGGACTTCATAATAAATGG--- +s Microtus_pennsylvanicus.CM073169.1 125811949 166 + 145900345 TGTTTCtgtcttgaaaattc-aaaaacaacaattaatgtgttt-aatatttatcaaaatagtgtaagtcaaatTTTatcttactattttaaaaactaaaaaa---aaATGG--CTCATTCCT-AT-CATTGGACTCTCC-AAATCTGTAGGCTTCAAAGGACTTCATAATAAATGG--- +s Neotoma_floridana.CM113784.1 161784805 164 + 180930968 TGTTTCTGTCTTGAAATTTC-AAACACAGGAATCAATGTTTTT-AATATTTACCAAAATAGTTTAAGTCAAATTT-ATCTTATTACTTTCAAAACTAAGAAA---AA--TG-GCTCATTTCT-AT-AATAGGACTCTCAG-AATCTGTAGGGTTCAAAGGATTGAATAGTAAATGG--- +s Onychomys_torridus.LR877189.1 162645586 164 - 167895514 TATTTCTGTCTTGAAATTTC-AAAAACAGCAATCAATGTTTTC-GATAGTTACCAAAATAGCTGAAGTGAAATTT-ATCTTATTACTTTCAAAACTAATAAA---AA--TA-GCTCATTTCT-AT-AATAGGACACCCAG-AATCTGTAGGATTCTAAAGACTGAACAGTAAATGG--- +s Peromyscus_californicus.VALE03000005.1 30491372 164 - 169319897 TATTTCTGTCTTGAAGTTTC-AAAAACAGCaatcaatgttttc-aatagttaccaaaatagtttaagtcaaatTT-atcttattactttaaaaactaagaaa---aa--tg-gcTCGTTTCT-AT-AATAGGACTCTCAG-AATCTGTAGCGTTCTAAGGAATGAACAGTAAATGG--- +s Peromyscus_maniculatus_North.CM113196.1 38376219 164 + 173669867 TATTCCTGCCTTGAAATTTC-AAAAGCAGCAATCAATGTTTTC-AATAGTTACCAAAATAGTTTAAGTCAAACTT-ATCTTATTACTTTAAAAACTAAGAAA---AA--TG-GTTCATTTCT-AT-AATAGGACACTCAG-AATCTGTAGGGTTCTAAGGATTGAACAGTAACTGG--- + +a +s Mus_musculus.CM000994.3 3602158 290 + 195154279 TGGTGCACGTCATGGTCACCTGGCTGGGCTTGTACATCTACCAGAAACTGGACATCTTTGCCCTTAAAAATCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAAACTTTCTCTGCAGAACAACACCATAGGCATCTATCAGCTGGCCAAGGCCATGACCACGCCGGTTATCATAGCCATCCAAACCTTCTGGTACCAGAAGAGATTCTCC------ATCAGAATACAGCTCACCCTGATTCCTATTACTGTA---------GG------TGTAATCCTAAAT +s Rattus_norvegicus.CM070397.1 81778891 276 - 137014596 TCGTGCACTTCGTGGTCACCTGGCTGGGCTTGTACGTCTGCCAGAAACTGGACATCTTTGCCCCCAAGAGTCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTTTCTCTGCAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACGACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAGAAACGATTCTCC------ATCAGAATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------T------------ +s Phodopus_campbelli.chr11 9440789 275 - 43870749 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCCCCCAAAAGTCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGCGGCTTTGTGGTCTTCACCAACCTTTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAGAAGAGCTTCTCCATCCGG------ATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------------------- +s Phodopus_sungorus_new.chr13 33378444 275 - 34507086 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCCCCCAAAAGTCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGCGGCTTTGTGGTCTTCACCAACCTTTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAGAAGAGCTTCTCCGTCCGG------ATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------------------- +s Phodopus_roborovskii.JAEOAF010000043.1 9199349 275 - 14393737 TGGTGCACTTCGTGATCACCTGGCTGGGCTTGTACATCTGCCAGAAACTGGACATCTTCGCCCCCAAAAGTCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGCGGCTTTGTGGTCTTCACCAACCTTTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAGAAGAGCTTCTCCATCCGA------ATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------------------- +s Cricetus_cricetus.OZ199630.1 147280374 275 + 245507489 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAAATGGACATCTTTGCCCCAAAAAGTCTGCCGCTCTCCAAGATCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTTTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCAGTGATCATAGCCATCCAGACCTTCTGGTACCAGAAGAGCTTCTCCATCCGA------ATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------------------- +s Cricetulus_griseus.NW_023276806.1 138612282 273 + 275698159 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAAATGGACATCTTTGCCCCCAAAAGTCTGCCGCTCTCCAAGATCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTTTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACTACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAGAAGAGCTTCTCCATCCGA------ATACAGCTCACGCTGGTGAGTA--GCTGGA------------------------------ +s Dicrostonyx_torquatus.JAQHUG010000343.1 1347017 269 - 2031119 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCTCCCAAAAGTCTGCCTCTCTCCAAGCTCTTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTTTCTCTGCAAAACAACACTATAGGTACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATACAGCTCACAATGGTGAGT--------A---------GC------------------- +s Arvicola_amphibius.NC_052063.2 42096221 269 + 42645200 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCTCCCAAAAGTCTGCCTCTCTCCAAGCTGCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTTTCTCTGCAAAACAACACTATCGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATACAGCTCACGATGGTGAGT--------A---------GC------------------- +s Microtus_pennsylvanicus.CM073186.1 3840003 269 + 39339920 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCTCCCAAAAGTCTGCCTCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTTTCTCTGCAAAACAACACTATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATACAGCTCACGATGGTGAGT--------A---------GC------------------- +s Microtus_arvalis.JAHRII010003550.1 36531763 269 + 36923249 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCTCCCAAAAGTCTGCCtctctccaagctcctcctcctggccctcagcttctgtggctttgtggtcttcaccaaccttTCTCTGCAAAACAACACTATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATACAGCTCACGATGGTGAGT--------A---------GC------------------- +s Neotoma_floridana.CM113806.1 12038210 275 + 53001843 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCCCCCAAAAGTCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTTTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCAGTGATCATAGCCATCCAGACCTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------------------- +s Onychomys_torridus.LR877207.1 12064208 275 - 48895294 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCCCCGAAAAGTCTGccgctctccaagctcctcctcctggccctcagcttctgtggctttgtggtcttcACCAACCTCTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------------------- +s Peromyscus_californicus.VALE03000023.1 11592836 275 - 46423440 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCCCCCAAAAGTCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTCTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACCTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATACAGCTCACGCTGGTGAGTA--GCTGCA---------GC------------------- +s Peromyscus_maniculatus_North.CM113212.1 52752033 284 + 54680123 TGGTGCACTTCGTGGTCACCTGGCTGGGCTTATACATCTGCCAGAAACTGGACATCTTTGCCCCCAAAAGTCTGCCGCTCTCCAAGCTCCTCCTCCTGGCCCTCAGCTTCTGTGGCTTTGTGGTCTTCACCAACCTCTCTCTACAGAACAACACCATAGGCACCTATCAGCTGGCCAAGGCCATGACCACGCCGGTGATCATAGCCATCCAGACTTTCTGGTACCAAAAGAGCTTCTCCGTCAGA------ATacagctcacgctggtgagta--gctgcaGCTGCTGCTgc------------------- + +a +s Mus_musculus.CM000994.3 3602448 324 + 195154279 TCTTATTACGACCTGAAGTTTCATTCCCTTGGAATGGTGTTTGCCACCCTTGGTGTGGTAGTCACGTCCCTTTAACAAGTGTGGGTCGGAGCCAAGCAGCACGAACTGCAGGTGAACTCCATACAGCTACTTTACTACCAGGCCCCCGTGTCTTCCACCATGCTGCTGGTGGCTGTGCACTTCTTCGAGCCAGTGTTTGCAGAGGGAGTAGTATTTGGTTCCTGGTCAGTTTCTGCTTTGCTTATGGTACTGCTGTCTGGAGTAATCGCTTTCATGGTGAACTTATCGATTTACTGGATCATTGGG---------AACATGT--CACCAGTCACC +s Rattus_norvegicus.CM070397.1 81784041 214 - 137014596 CGTCA---------------------------------------------------------------------------GTGGGTCGGAGCCAAGCAGCATGAACTGCAGGTGAACTCCATGCAGCTGCTTTACTACCAGGCCCCCATGTCTTCGGCCATGCTGCTGGTGGCTGTGCCCTGCTTTGAGCCAGTGTTTGCAGAGGGAGGAATATTTGGTCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTTCTAATTGTCCTGGTATCTAATTGAACTGGTACAGTCATCTTT +s Phodopus_campbelli.chr11 9445621 210 - 43870749 TTTCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCATGAGCTCCAAGTGAATTCCATGCAGCTGCTGTACTACCAGGCACCTATGTCCTCTGCCATGCTGCTGGTGGCTGTGCCCTTCTTTGAGCCAGTGTTCGGAGAGGGAGGAATATTTGGTCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTCCTGGTTGTTTTGGTGCCTAA--GAAACTA--TAGTAGCCACT +s Phodopus_sungorus_new.chr13 33383094 210 - 34507086 TTTCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCATGAGCTCCAAGTGAATTCCATGCAGCTGCTGTACTACCAGGCACCTATGTCCTCTGCCATGCTGCTGGTGGCTGTGCCCTTCTTTGAGCCAGTGTTCGGAGAGGGAGGAATATTTGGTCCCTGGTCAGTGTCTGCTTTGGT----------------------------------------------AAGTCCTGGTTGTTTTGGTGCCTAA--GAAACTA--TAGTAGCCACT +s Phodopus_roborovskii.JAEOAF010000043.1 9203660 210 - 14393737 TTTCA---------------------------------------------------------------------------GTGGGTAGGAGCCAAGCAGCATGAGCTCCAAGTGAATTCCATGCAGCTGCTGTACTACCAGGCGCCCATGTCGTCTGCCATGCTGCTGGTGGCAGTGCCCTTCTGTGAGCCAGTGTTTGGAGACGGAGGAATATTTGGTCCCTGGTCGGTCTCTGCTTTGGT----------------------------------------------AAGTCCTGGTTGTTTTGGTTTCTAA--GAGACGA--TGGCAGTCACT +s Cricetus_cricetus.OZ199630.1 147285152 210 + 245507489 TTTCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCATGAACTACAAGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCCTCTGCCATGCTGCTGGTGGCTGTGCCCTTCTTTGAGCCCGTGTTTGGAGAGGGAGGAATATTTGGCCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTCCTGATTGTTTTGGTGCCTAA--GAAACTG--CAGTACTCACA +s Cricetulus_griseus.NW_023276806.1 138617169 210 + 275698159 TTTCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCATGAACTACAAGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCCTCTGCCATGCTGCTGGTGGCTGTGCCCTTCTTTGAGCCAGTGTTTGGAGAGGGAGGAATATTTGGCCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTCCTGATTGTTTTGGTGCCTAA--GAAACTG--CAGTAGGCACA +s Dicrostonyx_torquatus.JAQHUG010000343.1 1352324 210 - 2031119 TTTCA---------------------------------------------------------------------------GTGGGTGGGCGCCAAGCAGCATGAATTGCAGGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCTTCCGCCATGCTGCTGGTGGCTGTGCCTTTCTTTGAGCCAGTGTTCGGAGAGAGAGGAATATTTGGTCCCTGGTCAGTTTCTGCTCTGGT----------------------------------------------AAGTCCTGATTGTTTTGGTGCCTAA--GAAATTG--CAGAAGTCATC +s Arvicola_amphibius.NC_052063.2 42101110 210 + 42645200 TTTCA---------------------------------------------------------------------------GTGGGTGGGCGCCAAACAGCATGAATTGCAGGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCTTCCGCCATGCTGCTGGTGGCTGTGCCTTTCTTTGAGCCAGTGTTTGGAGAGGGAGGAATATTTGGTCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTCCTGATTATTTTGGTGCCTAA--GAAATTG--CAGAAGTCATT +s Microtus_pennsylvanicus.CM073186.1 3845206 210 + 39339920 TTTCA---------------------------------------------------------------------------GTGGGTGGGCGCCAAGCAGCATGAACTGCAGGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCTTCCGCCATGCTGCTCCTGGCTGTGCCTTTCTTTGAGCCCGTGTTTGGAGAGGGAGGAATATTTGGTCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTCCTGATTGTTTTGGTGCCTAA--GAAACTG--CAGAAGTCATT +s Microtus_arvalis.JAHRII010003550.1 36537046 210 + 36923249 TTTCA---------------------------------------------------------------------------GTGGGTGGGCGCCAAGCAGCATGAGCTGCAGGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCTTCCGCCATGCTACTGGTGGCTGTGCCTTTCTTTGAGCCCGTGTTTGGAGAAGGAGGAATATTTGGTCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTCCTGATTGTTTTGGTGCCTAA--GAAACTG--CAGAAGTCATT +s Neotoma_floridana.CM113806.1 12043716 210 + 53001843 TTGCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCATGAATTGCAAGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCTTCTGCCATGCTGCTGGTGACTGTGCCCTTCTTTGAGCCAGTGTTTGGAGAGGGAGGAATATTTGGTCCCTGGTCAGTTTCTGCTTTGGT----------------------------------------------AAGTCCTGATTATTTTGGTGCCTAA--GAAATTG--CAGTAGCCACT +s Onychomys_torridus.LR877207.1 12070121 210 - 48895294 TTTCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCATGAACTGCAGGTGAATTCCATGCAGCTGCTCTACTACCAGGCTCCCATGTCCTCCGCCATGCTGCTGGTGGCCGTGCCCTTCTTTGAACCAGTGTTCGGAGAGGGAGGAATATTCGGACCCTGGTCAGTTTCTGCTCTGGT----------------------------------------------AAGTCCTGACGGTTTCGGTGCCTAA--GAAGTTG--CGGTTGACACT +s Peromyscus_californicus.VALE03000023.1 11598801 210 - 46423440 TTTCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCACGAACTGCAAGTGAATTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCCTCCGCCATGCTGCTGGTGGCCGTGCCCTTCTTTGAGCCCGTGTTCGGAGAGGGAGGAATATTTGGTCCCTGGTCGGTTTCTGCTCTGGT----------------------------------------------AAGTCCTGATTGCTTTGGTGCCTAA--GGAATTG--CAGTGGCCACT +s Peromyscus_maniculatus_North.CM113212.1 52757950 210 + 54680123 TTTCA---------------------------------------------------------------------------GTGGGTGGGAGCCAAGCAGCACGAACTGCAAGTGAACTCCATGCAGCTGCTGTACTACCAGGCTCCCATGTCCTCCGCCATGCTGCTGGTGGCCGTGCCCTTCTTTGAGCCGGTGTTCGGAGAGGGAGGAATATTTGGCCCCTGGTCAGTTTCTGCTCTGGT----------------------------------------------AAGTCCTGGTTGTTTGGGTGCCTAA--GAAACGG--CAGTAGCCACT + +a +s Mus_musculus.CM000994.3 3602772 3 + 195154279 TAC +s Rattus_norvegicus.CM070397.1 81784550 3 - 137014596 gga +s Phodopus_campbelli.chr11 9445927 1 - 43870749 --C +s Phodopus_sungorus_new.chr13 33383400 1 - 34507086 --C +s Phodopus_roborovskii.JAEOAF010000043.1 9204035 1 - 14393737 --c +s Cricetus_cricetus.OZ199630.1 147285526 1 + 245507489 --c +s Cricetulus_griseus.NW_023276806.1 138617544 1 + 275698159 --C +s Neotoma_floridana.CM113806.1 12044092 1 + 53001843 --t +s Onychomys_torridus.LR877207.1 12070497 1 - 48895294 --T +s Peromyscus_californicus.VALE03000023.1 11599177 1 - 46423440 --T +s Peromyscus_maniculatus_North.CM113212.1 52758327 1 + 54680123 --T + +a +s Mus_musculus.CM000994.3 3602775 681 + 195154279 AACATGTTCGGACACTTCAAGTTCTGCATCACTCTCTGCGGAGGATACATTCTATTTAAGGATCCCCTGTCGGTAAACCAGGGACTCCGCATCCTATGCACATTCTTCGGTATCCTCACCTACACCCACTTTAAACTCAGTGATCAGGAAGACAGTAAGAGAGAGTAAGCTGGTCCAACGGCCCT-AGCTTTGGCAGAGA-AGAGAAAGGCAAATCATGAGGATGAG-------AAGTGACAAACATGTA---CGTCTATGGAAAAGCTCCC-TGTAGAATTCAGTCAG-TGACGTACAAGAAGTATTTTCTTCGCAGATAATGACAATTTTTAGTACCTTTTATCGATTTAA----TTTTT------AAAAA--CCAAATGGT----TCAGGCTTCTC---AAAGACTTA-ATGGACctgg-agctggcactcagtagtagagtgtttgcctcacacgcg--------tgagaccctggcctcaaactccaacaccGTCAAAGAAATGCTAGCAACTGCCTAGAAGATGGCTT-CATCCTGCCAGCTGAAGCACGGGAGAGTTCGGCTCCATAGGAATGACA-------GAGCGCTCCgtcgctgcagctgtgctgtggct-gtcgg-----tgtccagtgtccctgtTCTTCCAGGGACGGGCACGCCAG---CCT---GCTGGCTCACCCAGGCCTTCCTGTCAAGAAGTTGAAGTTCTTTGTCGTTCGGAC------------AGGG--TGGTG +s Rattus_norvegicus.CM070397.1 81789753 629 - 137014596 AACATGTTCGGACATTTCAAGTTCTGCATCACCCTGTGTGGAGGCTACATCTTATTTAAAGATCCACTGTCGGTAAACCAAGGGCTCGGCATCCTATGCACACTGTTTGGTATCCTCGCCTATACCCACTTTAAACTCAGCGAACAGGAAGGCAGTA----AGAGTAAGCTGGTCCAGCGGCCCTGAGCTTTGACAGAGA-AGAGAAAATGGCACCAAAAGAAGGAG--------------AAACA--GC---AGTGCGTGGAAAAGCTCCC-TGCAGAATTCAGTCAT-TGATGTACAAGAAGTATTTTCTTCACAGATAATGA--GTTTTTAGTACCTTTTATAGATTTAATTTTTTTTT------TAAAC--CAAAATGGT----TCAGGCTTATCTTTAAAGGCttt-atgaacctgg-agctgacactcagtggtcgagtgtttgcctcacacgca--------tgagaccctggcctcaaaccccagcaccgtcaaggaaatgctc-------------------------------------aagtacaggaaagcTCAGTTCCCCTGGAATGACG------TGGACGCTCCAGTGCTGCAGCTGTGCTGTGGCA-GTCGG-----TGTCCG-------TGACCTCCCGGGGAAGGGCACGCCAA---ACTACCGCTGCCTCACCCATGCCTTCCCATCGG------GAAGTTCTTCGTCCTTCAGAC------------Aggg--tagtg +s Phodopus_campbelli.chr11 9452104 622 - 43870749 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGTGGAGGATACATCTTGTTTAAGGATCCACTGTCCGTTAACCAGGCACTCGGCATCTTGTGCACACTGTTCGGCATCCTCACCTATACCCACCTGAAACTCAGCGAACAGGAAGGGAGCA----AGAGTAAGCTGGTCCAGCGGCCCTGAGCTTGGCCGG-GGAAAAGAAAGCTGCCCCAGAAAGACTAGAGACGTCACATGTC-------C---------GTGGAAAAGCTAC------------------------------AAGTGTTTTCTTCTCACA-GATGAC---------TGCCTTCCATCAATTGGG----TTTGAT-----GGACCAACAAAATGGTCCAC----ACTTCTTTTCCAAGTCTTC-ATGAAGCTGC-AGCTCTCActcattgatggagtgtttgcctagcatgca--------tcagaccctgggctcaaaccccagCACCTTCAAAGAAATGATGTGAACTACATATAAGATAGTTT-CATCCTACCGACTAAAATATAGGGAAGTTCGGTCCCGTAGGAATAACAAA-G--------------CGCTACGGGCATGCCCTGGCTTGTCCGGTTGGCGGTCA-------CACTCTCCCCAGGAAGGGCACGTAGGACTTCT---GCCCTCTCA----------CCTACTGTGAGGTTTG--TTACATGGTAGTCCGTG------------TGTG--CAGTG +s Phodopus_sungorus_new.chr13 33389685 632 - 34507086 AACATGTTCGGGCACTTCAAGTTCTGCATCACCCTGTGTGGAGGATACATCTTGTTTAAGGACCCACTGTCCGTTAACCAGGCACTCGGCATCTTGTGCACACTGTTCGGCATCCTCACCTATACCCACCTGAAACTCAGCGAACAGGAAGGGAGCA----AGAGTAAGCTGGTCCAGCGGCCCTGAGCTTGGCCGG-GGAAAAGAAAGCTGCCCCAGAAAGACTAGAGACGTCACATGTC-------T---------GTGGAAAAGCTAC------------------------------AAGTGTTTTCTTCTCACA-GATGAC---------TGCCTTCCATCAATTGGG----TTTGCT-----AAAAC--CAAAATGGTCCAC----ACTTCTTTTCCAAGTCTTC-ATGAAGCTGC-AGCTCTCACTCAttgatggagtgtttgcctagcatgca--------tcagaccctgggctcaaaccccagcATCTTCAAAGAAATGATGTGAACTACATATAAGATAGTTT-CATCCTACCGACTAAAATATAGGGAAGTTCAG--CCGTAGGAATAACAAA-GGGTGTGCACTCCAGCGCTACGGGCGTGCCCTGGCTTGTCCAGTCGGCGGTCA-------CACTCTCCCCAGGAAGGGCACATAGGACTTCT---GCCCTCTCA----------CCTACTGTGAGGTTTG--TTATATGGTAGTCCGTG------------TGTG--CAGTG +s Phodopus_roborovskii.JAEOAF010000043.1 9208366 620 - 14393737 AACATGTTCGGACACTTGAAGTTCTGCATCACCCTGTGCGGAGGATACATCTTGTTTAAGGACCCACTGTCCGTTAGCCAGGGCCTCGGCATCTTGTGTACACTGTTCGGCATCCTCACCTACACCCACCTGAAACTCAGCGAACAGGAAGGGAGCA----AGAGTAAGCTGGTCCAGCGGCCCTGAGCACTCACGG-AA-A--------------AGAAAGACTAGGGACCTCAAATGCC-------TAAGTC---CGTGGAAAGGCTAC------------------------------AGGTGTTTTCTTCTCAAA-GATGAC---------TGCCTCCCATCAATTGGA----TTTGCT-----AAAAC--CAAAGTGGTCCAC----GCTTCTTTTCCAAGTCTTC-ATGAAgctgc-agctgtcactc--------agtgtttgccaagcatgca--------tcagaccctgggctcaaaccccagcaccttcaaagaaaTGATAGGAGCTACATGTAAGATGGTTT-CATCCTACTGACTAAAGTACAGGGAAGTTCGGTTCCCTAGGAATGGCAAA-GCGTGTGCGCTCCAGCGCTACAGGCGCGCCCTGGC--ATCCGGTCG-CGTTCA-------CGCTCTCCCCAGGAAGGGCACGTAGGACTTCT---GCCCTCTCATGCATGGTTTCCTACTGTGAGGGTTA--CTCTATGTTA----CTG------------TGta--cagtg +s Cricetus_cricetus.OZ199630.1 147291859 689 + 245507489 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGTGGAGGATACATTTTGTTTAAGGATCCACTGTCAGTTAACCAGGGACTTGGCATCTTATGCACACTGTTTGGCATCCTCACCTATACCCACTTTAAACTCAGCGAACAGGAAGGAAGCA----AGAGTAAGCTGGTTCAACGGCCCTGAGCTCTACCAG-GG-AAAGAAAGCTTCACCAAAAGGATGAGAGACGTCAAATGTC-------TAAGTC---CATGGAAAAGCTACCATGCAAAGTTCACT-GGATGACTTACCAGAAGTATTTTCTTCACAAA-GATGACTGCTTTTAATACCTTTCATCGATTTGA----TTTTCT-----AAAAC--CAAAATGGTCCAC----GCTTCCTTTTAAAGTCTTT-Atgaagcggc-agctgtcactcattggtagagtgtttgcctagcatgcaCTAGCATAtcagactctgggctcagaccccagcactttcaaagaaattatatgaactgcatataaaatggttt-cattctaccaactaaaaTATAGGAAAGTTTGGTTCCATAGGAATGACGAA-GGGTGCACACTCCAATGCTGCAGGTGTGCCCTGGCTTGTCGG-----TGGTCA-------CACTCTCCCCAGGAAGGGCACATAGGACTTCT---GTCGACTCATCCATGGTTTCCTACTATGAGGTTTA--TTATAtgttagt-ggtc------------tgtg--tggtg +s Cricetulus_griseus.NW_023276806.1 138624061 690 + 275698159 AACATGTTCGGACATTTCAAGTTCTGCATCACCCTGTGCGGAGGATACATTCTGTTTAAGGATCCACTGTCAGTTAACCAGGGGCTTGGCATCTTATGCACACTGTTCGGCATCCTCACCTATACCCACTTTAAACTCAGCGAACAGGAGGGAAGCA----AGAGTAAGCTGGTCCAACGGCCCTGAGCTCTACCAG-AG-AAAGAAAGCTTCACCAAAACGATGAGAGACGTCAAATGTC-------TACGTC---CATGGAAAAGCTATCATGCAAAGTTCACT-GGATGACTTACCAGAAGTATTTTCTTCACAAA-GATGACTGCTTTTAATACCTTTCATCGATTTGA----TTTTCT-----AAAAC--CAAAATGGTCCAC----GCTTCCTTTTAAAGTCTTT-ATGAAGCggc-agctgtcactcattgatagagtgtttgcctagcatgcaCTAGCATAtcagaccctgggctcagaccccagcactttcaaagaaatgatatgaactgcatataagatggttt-cattctaccaactaaaaTATAGGAAAGTTCGGTTCCATAGGAATGACGAA-GGGTGGGCACTCCAATGCTGCAGGTGTGCCCTGGCTTGGCGG-----TGGTCA-------CACTCTCCCCAGGAAGGGCACATAGGACTTCT---GCCTACTCATCCATGGTTTCCTACTGTGAGGttta--ttatatgttagtcagtc------------tgtg--tggtg +s Dicrostonyx_torquatus.JAQHUG010000254.1 1866509 647 - 2996520 GACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGCGGAGGATACCTTTTATTTAAGGACCCACTGTCGGTTAACCAAGGACTCGGCATCTTATGTACGCTGTTTGGCATCCTCACCTACACCCACTTTAAACTCAGCGAGCAGGAAGGAACTA----AGAGTAAGCTGGCCCAACGGCCCTGAGCTTTggcggaga-aaagaaggctgtaccaaagagatgag-------aa-tgtg-------gaAGTc-----tggaaaagctaccAtgcagagttctccag--tgacttacaagaagtattttcttca-aaa-tatgacgactgttggtacttttcatcagtttaa----tgttcT-----aaaac--caaa-tggtTAGC----gcttctcTTGaaagacttt-ctgaagct-g-agctgtcactcagtggcgcagcattggcctcgcatgca--------taagaccctgggctcaaactccggcgcctacaaagaaatgaCACGAACTGCAT----------------TCTGACAACTGAAATACAGAACAGTTCGGTTCTGTAGGAATGACAGA-GGGTGCACAGTCCAGTGCTGCAGGCGTGCCCTGGCTTGTTTG-----TGGTAA-------CACGCTCCCCAGGAAGGGCACATAGGACTcct---gcctgctaa-ccatggtttcccactgtgaggttta--gtgtatgttagtcac--------------tgtg--t---- +s Arvicola_amphibius.NC_052063.2 42109121 613 + 42645200 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGTGGAGGATACATTTTATTTAAGGACCCGCTGTCAGTTAACCAGGGACTTGGCATTTTATGTACGCTGTTTGGCATCCTCGCCTATACCCACTTTAAACTCAGCGAACAGGAAGGAAGTA----AGAGTAAACTGGTCCAACGGCCCTGAGCTTTGGCGGAGA-AAAGAAGGCTGAG--AAAAAGATGAG-------AAATGTG-------TAAGTC---TGTGGAAAAGCTACCATGCAACATTCCCCCAG-TGACTTACAAGAAGTATTTTCTTCA-AAA-TGTGACGACTGTTAGTACTTTTCATCGATTTAA----TTTTCT-----AAAAC--CAAAATGGTTAGC----GCTTCTCTTCAAAGACTTT-CTG------G-AGCTGCCACTCAGCAGTGCAGCTCC-----------------------AGACCCTGGGCTCAAACTCCGGCGCCTACGAAGAAGTGGCACGGACTGCGTGGAAGGTGG-TTTCATCCTGCCAGCTGAGACACAGGAAAGTTGTGTTCTGTGGGAAGGACAAA-GGGTGCACGCTCCAGTGCTGCAGGCA---------------G-----TGCT---------------------------TACATAGCACTCCT---GCCCGCTAA-C-AGGGTTTCCCACTGTGAGGTGTA--GTATACGTTAgtcac--------------tgtg--cagca +s Microtus_pennsylvanicus.CM073186.1 3854171 591 + 39339920 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGCGGAGGATACGTTTTATTTAAGGACCCGCTGTCGGTTAACCAGGGACTCGGCATCTTATGTACGCTGTTCGGCATCCTCACCTATACCCACCTTAAACTCAGCGAACAGGAAGGAAGTA----AGAGTAAACTGGTCCAACGGCCCTGAGCTTTGGCGGAGA-AAAGAAGGCTGTACCAAAAAGACGAG-------AAGTGTG-------TAACTC---AGTGGAAAAGCAGCCGTGCAACGTTCCCCCAG-TGACTTACA--------TTTCTTCA-AAA-TGTGACGACTGGTAGTACTTTTCATCGATTTAA----TTTTCT-----AAAAC--CAAAATGGTTAGC----GCTTCTCTTCAAAGACTTTTCTGAAgct-gGaactgtcactcagcggtgcagcgtttgcctcgcatgca--------gaagaccctgggctcaaaccccggcgcctacaaagaa-TGACACAGACTGCGTGTAAGATGG-TT-CATTCTACCAACTGAAGTACAGGAAAGTTGGGTTCCT-----GTG-----------------------CTGCAGACATGCCCTGGCTCCTCTG-----TGCT---------------------------TACAGAGGACTCCT---GCCCGCTAA-C-ATGGTTTCCCACTGTG----------------------------------------TG--CAGCT +s Microtus_arvalis.JAHRII010003550.1 36542465 642 + 36923249 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGCGGAGGATACGTTTTATTTAAGGACCCGCTGTCAGTTAACCAGGGACTCGGCATCTTATGTACACTGTTTGGCATCCTCGCCTATACCCACTTTAAACTCAGCGAACAGGAAGGAAGTA----AGAGTAAGCTGGTCCAACGGCCCTGAGCTTTGGCGGAGA-AAAGAAGGCTATACCAAAAAGATGGG-------AAGTGTG-------CAAGTC---CGTGGAAAAG-AACCGTGCAGCGTTCCC--AG-TGAATTACAGGAGG-ATTGTCTTCA-AAA-TGTGACGACTGCTAGTACTTTTCATCTATTTAA----TATTCT-----AAAAC--CAAAATGGTTAGC----GCTTCTCTTCAAAGACTTT-CTgaagct-gGagctccaactcagcggtgcagcgtttgcctcgtgtgca--------gaagaccctgggctcaaactccggcgcctacaaagaa-tgacacagactgcgtgtaagacag-tt-cATTCTACCAACTGAAGTACAGGAAAGCTGGGTTCCT-----GTGACAAA-GGGTGCACACTCCAGTGCTGCAGGCGTGCCCTGGCTCATCTG-----TGCT---------------------------TACATAGGACTCCT---GCCCGCTAA-C-ATGGTTTCCCACTGTGAGGttta--gtatatgttagtcac--------------tgtgTGcagct +s Neotoma_floridana.CM113806.1 12052782 681 + 53001843 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGCGGAGGCTACATTTTGTTTAAGGATCCACTGTCGTTTAACCAGGGACTCGGCATCTTATGCACACTGTTCGGCATCCTCACCTATACCCACTTTAAACTCAGCGAACAGGAAGGAAGTA----AGAGTAAGCTGGCCCAACGGCCCTGAGCTTTGGCGGAGA-AAAGAAAGCTGCACCCAAAAGATGAG-------AAATGTCAAATGTATAAGTC---CCTGGAGAAGCCATCATGCAAAGTTCACTCAG-TGACTTACAAGAAGTATTTTCTTCACAAG-TACGACTGTTTTTAGTACCTTTCATCGATTTAA----TTTTCT-----AAAAC--CAAAATGATTTGT----GCTTCTCTTCAAAGACTTT-ATgaagctgg-agctgtcactcaatggtagagtgtttgcctagcatgca--------taagaccctgggctcaaagtccagcaccttaaaaggagtgatacaaactacatacaagatggttt-caTTCTACCAACCAAAATGTAGAAAAGTTTGGTTCCATAGGAATGACAAAGGGGTGCACACTTCAGTGCCGCAGGTGTGCCCCGGTTTG---G-----TGGTTA-------TGCTCTCCCCAGGAAGGGCACATAGGACTCCT---GCCCACTCATCCATGGTTTCCCACTGTGAGGTTAA--TTATATGTTGGTCAGTG------------TAtg--cagtg +s Onychomys_torridus.LR877207.1 12077813 678 - 48895294 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGCGGAGGATACATTTTATTTAAGGATCCACTGTCCTTTAACCAGGGACTCGGCATCTTATGCACACTGTTCGGCATCCTCACCTATACCCACTTTAAACTCAGCGAACAGGAAGGAAGTA----AGAGTAAGCTGGTCCAGCGGCCCTGAGCTTTGGCGGAGA-AGAGAAAGCCGCACCAAAAAGATGGG-------AAATGC-------CTAAATC---ACTGGAGAGGCCACCATGCAGAGTTCACTCAG-TGACTTACAAGAAGTATTGTCTTCACAAA-TCTGAATGTTTTCAGTACCTTTCATCGACGTAA----TTTTCTAAATTAAAAC--CAAAATGGTTGGA----GCTTCTCTTCGAAGACTTT-ATGAAGCagg-agccgtcactcagcggtagagtgtctgcccagcatgtg--------taaggcccca-gctcagactccgacaccTTCAAAGAGACAATACATACTCCATATGAGATGGTTT-CATTCTGCCAACCAAAATATAGGAAGGTTCGGGGCCATAGGAACAACAAAGGGCTGCATACTCCAGTGCCGCAGGTATGCCCCGGCTTG---G-----GGGTCA-------CGCTCTCCCCAGGAAGGGCACATAGGACTCCT---GCCCACTCACCCCCGGTTTCCCACTGTGAGGTTTA--TTATATGCTAGTCAGTG------------TAtg--cagtg +s Peromyscus_californicus.VALE03000023.1 11607620 678 - 46423440 AACATGTTCGGACACTTCAAGTTCTGCATCACCCTGTGCGGAGGATACATTTTATTTAAGGATCCGCTGTCGTTTAACCAGGGACTTGGCATCTTATGCACACTGTTCGGCATCCTCACCTATACCCACTTTAAACTCAGCGAACAGGAAGGAAGTA----AGAGTAAGCTGGTCCAGCGGCCCTGAGCTTTGACGGAGA-AAAGAACGCTG-ACCAAAAAGCTGGG-------AAATGG-------ATAAGTC---CCTGGAGAAGCCACCATGCAGAGTTCACTCAG-TGACTTACGAGAAGTATTTTCTTCACAAA-T----CTGTTTTTAGTACCTTTCATCGACTTAA----TTTTCTAAA-----AC--CAAAATGGTTCGA----GCTTCTCTTCAAAGACTGC-ATGAAgctgg-agctgtcactcagcggtagagtgtctgctccgcatgca--------tcaggccccgggctccaactccagcaccttcaaagaaaCGATACATGCTACATAGGAGATGGTTT-CATTCTGCCAGCCAAAATGTAGGAAGGTTCGGGTCCATAGGAATGACAAAGGGCTGC-CACTC-AGTGCTGCAGGTGTG-CCCGGCTTG---G-----TGGCCA-------CGCTGTCCCCAGGAAGGGCACATAGGACTCCT---GCCCACTCATCCCTGGTTTCCCCCTGTGAGGTTTA--TtctatgttagtcagtgTTTATTCTATGTtatg--ccgtg +s Peromyscus_maniculatus_North.CM113212.1 52764572 667 + 54680123 AACATGTTTGGACACTTCAAGTTCTGCATCACCCTGTGCGGAGGCTACATTTTATTCAAGGATCCACTGTCGTTTAACCAGGGACTCGGCATCTTATGCACGCTGTGTGGCATCCTCGCCTATACCCACTTTAAACTCAGCGAACAGGAAGGAAGTA----AGAGCAAGCTGGTTCAGCGGCCCTGAGCCTCGGCGGAGA-AAAGAACGCTGCACCAAAAAGCTGGG-------AAATGT-------ATAAGTC---CCTGGAGAGGCCACCATGCAGAGTCCACTCAG-TGACTTACGGGAAGTATTTTCTTCACAAA-TCTGACTGTTTTTAGTACCTTTCATCGACTTAA----TTTTCTAAA-----AC--CAAAATGGTTCGA----GCTTCTCTTCAAAGACTTT-ATgaagctag-agcagtcactcagtggtagagtgtctgcctagtatgca--------taaggccccgggctcaaactccagccccttcaaagaaaCGATACATACTACATCTGAGATGGTTT-CATTCTACCAACCAAAATATAGGAAGGTCCGGGTCCATAGGAATGACAAAGGGCTGC-CACTC-AGTGCTGCAGGCGTGCCCCGGCTTG---G-----TGGCCA-------CGTTGTCCCCAGGAAGGGCACACAGGACTCCT---GCCCACTCATCCCTGG-TTCCCCCTGTGAGGttta--ttatgtgtt----agtg------------tatg--cagtg + diff --git a/tests/loci-classes.maf.block.idx b/tests/loci-classes.maf.block.idx new file mode 100644 index 0000000..d1600fa --- /dev/null +++ b/tests/loci-classes.maf.block.idx @@ -0,0 +1,7 @@ +# mafutils-index format=2 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 +CM000994.3 3602448 324 335 385 15 10577 16497 +CM000994.3 3602772 3 3 51 11 16497 17157 +CM000994.3 3602775 681 758 808 15 17157 29422 diff --git a/tests/loci-classes.maf.scaffold.idx b/tests/loci-classes.maf.scaffold.idx new file mode 100644 index 0000000..3518fd5 --- /dev/null +++ b/tests/loci-classes.maf.scaffold.idx @@ -0,0 +1,2 @@ +# mafutils-index format=2 maf=loci-classes.maf compression=none size=29422 mtime=1789746613.0 hash=md5:e2a9c32b01c9e0b4f798534f8815d6a2 +CM000994.3 16 29422 diff --git a/tests/test_loci.py b/tests/test_loci.py new file mode 100644 index 0000000..b8df17a --- /dev/null +++ b/tests/test_loci.py @@ -0,0 +1,613 @@ +import csv +import os +import subprocess +import sys +from collections import defaultdict + +import pytest + +from mafutils.fetch import (ELEMENT_LOCI_HEADERS, chunkGaps, classifyChunks, + fastaScaffoldField, summarizeElementLoci) + +TEST_DIR = os.path.dirname(__file__) +REPO_ROOT = os.path.abspath(os.path.join(TEST_DIR, "..")) + +MAF_FILE = os.path.join(TEST_DIR, "example.maf") +INDEX_FILE = os.path.join(TEST_DIR, "example.maf.block.idx") +BED_FILE = os.path.join(TEST_DIR, "example.bed") + +# 6 real blocks (15 species) covering two real regions chosen because they +# exhibit the rare classes that tests/example.maf cannot produce: +# multi_strand and multi_scaffold. Extracted read-only from production data. +CLASSES_MAF = os.path.join(TEST_DIR, "loci-classes.maf") +CLASSES_INDEX = CLASSES_MAF + ".block.idx" +CLASSES_REGIONS = [ + ("CM000994.3", 3541891, 3541950, "multistrand"), + ("CM000994.3", 3602239, 3602958, "multiscaffold"), +] + + +def run_fetch(args): + return subprocess.run( + [sys.executable, "-m", "mafutils", "fetch"] + args, + check=False, cwd=REPO_ROOT, capture_output=True, text=True, + ) + + +def write_bed(path, rows): + with open(path, "w") as fp: + for r in rows: + fp.write("\t".join(str(x) for x in r) + "\n") + + +def read_loci(path): + with open(path) as fp: + return list(csv.DictReader(fp, delimiter="\t")) + + +def headers(fa_path): + with open(fa_path) as fp: + return [ln.strip() for ln in fp if ln.startswith(">")] + + +# --------------------------------------------------------------------------- +# classifyChunks / fastaScaffoldField (unit level) +# --------------------------------------------------------------------------- + + +def chunk(scaffold="s1", start=0, size=10, strand="+", src_size=1000, bi=1): + """bi = block_index, the block ordinal this chunk came from.""" + return {"scaffold": scaffold, "start": start, "size": size, + "end": start + size, "strand": strand, "src_size": src_size, + "block_index": bi} + + +def test_single_chunk_is_single(): + assert classifyChunks([chunk()]) == ("single", 0) + + +def test_adjacent_chunks_are_contiguous(): + chunks = [chunk(start=0, size=10), chunk(start=10, size=5)] + assert classifyChunks(chunks) == ("contiguous", 0) + + +def test_gap_makes_it_split_and_reports_the_gap(): + chunks = [chunk(start=0, size=10), chunk(start=60, size=5)] + assert classifyChunks(chunks) == ("split", 50) + + +def test_largest_magnitude_gap_wins_and_sign_is_kept(): + """ + Overlapping chunks (negative gap) mean something different from a plain + gap, so the sign must survive rather than being abs()'d away. + """ + chunks = [chunk(start=0, size=10), chunk(start=5, size=10), chunk(start=100, size=5)] + cls, gap = classifyChunks(chunks) + assert cls == "split" + assert gap == 85 # 100 - (5+10), the largest magnitude + + overlapping = [chunk(start=0, size=50), chunk(start=10, size=5)] + cls, gap = classifyChunks(overlapping) + assert cls == "split" + assert gap == -40 # signed, not 40 + + +def test_differing_scaffolds_beat_other_signals(): + chunks = [chunk(scaffold="s1", start=0, size=10), chunk(scaffold="s2", start=10, size=5)] + assert classifyChunks(chunks) == ("multi_scaffold", None) + + +def test_differing_strands(): + chunks = [chunk(start=0, size=10, strand="+"), chunk(start=10, size=5, strand="-")] + assert classifyChunks(chunks) == ("multi_strand", None) + + +def test_scaffold_field_dedupes_preserving_block_order(): + chunks = [chunk(scaffold="b"), chunk(scaffold="a"), chunk(scaffold="b")] + assert fastaScaffoldField(chunks) == "b,a" + assert fastaScaffoldField([chunk(scaffold="")]) == "" + + +# --------------------------------------------------------------------------- +# The scaffold must be in the header regardless of dedupe / expected-species +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("extra", [ + [], + ["--fasta-dedupe", "most-seq"], + ["--expected-species", "human,chimp,gorilla"], + ["--fasta-dedupe", "most-seq", "--expected-species", "human,chimp,gorilla"], +]) +def test_scaffold_present_in_header_for_every_dedupe_mode(extra, tmp_path): + """ + The scaffold used to live only in the FASTA dict key, so keying by species + (which --fasta-dedupe and --expected-species both force) discarded it and + the header could not name a locus. + """ + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [("chr2", 5, 10, "r1")]) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([MAF_FILE, bed, "--index", INDEX_FILE, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id"] + extra) + assert result.returncode == 0, result.stderr + + hdrs = headers(os.path.join(out, "r1.fa")) + present = [h for h in hdrs if h.startswith((">human", ">chimp"))] + assert present, hdrs + for h in present: + assert h.startswith((">human.chr2", ">chimp.chr2")), h + + +def test_species_only_stays_bare(tmp_path): + """species-only exists so headers match tree tip labels; it never gains a scaffold.""" + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [("chr2", 5, 10, "r1")]) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([MAF_FILE, bed, "--index", INDEX_FILE, "-b", "id", "-o", out, + "-f", "-fh", "species-only", "--fasta-dedupe", "most-seq"]) + assert result.returncode == 0, result.stderr + assert sorted(headers(os.path.join(out, "r1.fa"))) == [">chimp", ">human"] + + +def test_species_coords_id_survives_a_bed_with_no_name_column(tmp_path): + """Regression: id_str was assigned conditionally but used unconditionally.""" + bed = os.path.join(str(tmp_path), "noid.bed") + write_bed(bed, [("chr2", 5, 10)]) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([MAF_FILE, bed, "--index", INDEX_FILE, "-o", out, + "-f", "-fh", "species-coords-id"]) + assert result.returncode == 0, result.stderr + # default --basename is 'coords', so the file is named scaffold-start-end + for h in headers(os.path.join(out, "chr2-5-10.fa")): + assert not h.endswith(" "), f"trailing space with no region id: {h!r}" + assert "id:" not in h + + +# --------------------------------------------------------------------------- +# The loci table +# --------------------------------------------------------------------------- + + +def test_loci_table_requires_fasta(tmp_path): + result = run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", + "-o", os.path.join(str(tmp_path), "out"), "--loci-table"]) + assert result.returncode != 0 + assert "--loci-table requires --fasta" in (result.stdout + result.stderr) + + +def test_loci_table_grain_and_class_consistency(tmp_path): + out = os.path.join(str(tmp_path), "out") + result = run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--loci-table"]) + assert result.returncode == 0, result.stderr + + rows = read_loci(os.path.join(out, "maf_fetch_loci.tsv")) + assert rows + + groups = defaultdict(list) + for r in rows: + groups[(r["region_scaffold"], r["region_start"], r["region_end"], r["species"])].append(r) + + for key, g in groups.items(): + # one row per contributing chunk, indexed 1..n, class constant per group + assert len({r["class"] for r in g}) == 1, key + assert len({r["max_gap"] for r in g}) == 1, key + assert int(g[0]["n_chunks"]) == len(g), key + assert sorted(int(r["chunk_index"]) for r in g) == list(range(1, len(g) + 1)), key + if len(g) == 1: + assert g[0]["class"] == "single", key + + +def test_loci_table_matches_the_examples_known_shape(tmp_path): + """ + 03-crossblocks-missing is a hand-checkable case: human's two chunks are 6bp + apart (split) while chimp's are adjacent (contiguous) in the same region. + """ + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [("chr1", 24, 33, "cross")]) + out = os.path.join(str(tmp_path), "out") + result = run_fetch([MAF_FILE, bed, "--index", INDEX_FILE, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--loci-table"]) + assert result.returncode == 0, result.stderr + + rows = read_loci(os.path.join(out, "maf_fetch_loci.tsv")) + by_sp = defaultdict(list) + for r in rows: + by_sp[r["species"]].append(r) + + assert by_sp["human"][0]["class"] == "split" + assert by_sp["human"][0]["max_gap"] == "6" + assert by_sp["chimp"][0]["class"] == "contiguous" + assert by_sp["chimp"][0]["max_gap"] == "0" + # src_size is emitted so forward coordinates are one subtraction away + assert all(int(r["src_size"]) > 0 for r in rows) + + +@pytest.mark.parametrize("processes", [1, 4]) +def test_loci_table_is_deterministic_across_process_counts(processes, tmp_path): + base = os.path.join(str(tmp_path), "p1") + assert run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-o", base, + "-f", "-fh", "species-coords-id", "--loci-table", "-p", "1"]).returncode == 0 + + out = os.path.join(str(tmp_path), f"p{processes}") + result = run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--loci-table", "-p", str(processes)]) + assert result.returncode == 0, result.stderr + with open(os.path.join(out, "maf_fetch_loci.tsv")) as a, open(os.path.join(base, "maf_fetch_loci.tsv")) as b: + assert a.read() == b.read() + + +def test_loci_table_does_not_change_the_fasta(tmp_path): + without = os.path.join(str(tmp_path), "without") + with_ = os.path.join(str(tmp_path), "with") + common = [MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-f", "-fh", "species-coords-id"] + assert run_fetch(common + ["-o", without]).returncode == 0 + assert run_fetch(common + ["-o", with_, "--loci-table"]).returncode == 0 + + for name in os.listdir(without): + if name.endswith(".fa"): + with open(os.path.join(without, name)) as a, open(os.path.join(with_, name)) as b: + assert a.read() == b.read(), name + + +# --------------------------------------------------------------------------- +# The rare classes, on real data +# --------------------------------------------------------------------------- + + +def test_real_fixture_produces_the_rare_classes(tmp_path): + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, CLASSES_REGIONS) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([CLASSES_MAF, bed, "--index", CLASSES_INDEX, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--fasta-dedupe", "most-seq", + "--loci-table"]) + assert result.returncode == 0, result.stderr + + rows = read_loci(os.path.join(out, "maf_fetch_loci.tsv")) + seen = {r["class"] for r in rows} + assert {"multi_scaffold", "multi_strand", "split", "contiguous"} <= seen, seen + + # max_gap is undefined across different source sequences + for r in rows: + if r["class"] in ("multi_scaffold", "multi_strand"): + assert r["max_gap"] == ".", r + + +def test_multi_scaffold_header_comma_joins_and_warns(tmp_path): + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [CLASSES_REGIONS[1]]) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([CLASSES_MAF, bed, "--index", CLASSES_INDEX, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--fasta-dedupe", "most-seq", + "--loci-table", "--verbose"]) + assert result.returncode == 0, result.stderr + + rows = read_loci(os.path.join(out, "maf_fetch_loci.tsv")) + multi = {r["species"] for r in rows if r["class"] == "multi_scaffold"} + assert multi, "fixture no longer exercises multi_scaffold" + + hdrs = headers(os.path.join(out, "multiscaffold.fa")) + for sp in multi: + h = next(x for x in hdrs if x.startswith(f">{sp}.")) + scaffolds = {r["src_scaffold"] for r in rows if r["species"] == sp} + assert "," in h, h + for scaf in scaffolds: + assert scaf in h, (scaf, h) + + assert "multi-scaffold" in (result.stdout + result.stderr) + + +# --------------------------------------------------------------------------- +# Zero-size chunks and the no_bases class +# --------------------------------------------------------------------------- + + +def test_all_zero_size_chunks_are_no_bases_not_single(): + """ + A species can have an s-line in a block whose aligned sequence falls + entirely outside the region once trimmed, leaving an all-gap row and a + zero-width header span. That used to classify as a clean `single` + (1,648 such rows on real data). + """ + assert classifyChunks([chunk(size=0)]) == ("no_bases", 0) + assert classifyChunks([chunk(start=0, size=0), chunk(start=50, size=0)]) == ("no_bases", 0) + + +def test_zero_size_chunks_do_not_create_phantom_gaps(): + """ + Real case (Dicrostonyx_torquatus, CM000994.3:10038193): a zero-size chunk + in one block plus 27 real bases in a later block. Measuring a gap from a + chunk containing no bases invented a 718bp `split`; the species really has + one aligned piece. + """ + chunks = [chunk(start=9949990, size=0), chunk(start=9950708, size=27)] + assert classifyChunks(chunks) == ("single", 0) + assert chunkGaps(chunks) == [] + + +def test_zero_size_chunk_is_still_emitted_as_a_row(tmp_path): + """Excluded from the arithmetic, but never hidden from the table.""" + from mafutils.fetch import writeLociRows + import io + buf = io.StringIO() + chunks = [chunk(start=100, size=0, bi=1), chunk(start=500, size=27, bi=2)] + writeLociRows(buf, {"scaffold": "chr1", "start": 0, "end": 1000, "id": "r1"}, + "sp", chunks, "single", 0) + lines = [l for l in buf.getvalue().splitlines() if l] + assert len(lines) == 2, lines + assert "\t0\t" in lines[0] # the zero-size chunk survives as a row + + +# --------------------------------------------------------------------------- +# Per-element summary columns +# --------------------------------------------------------------------------- + + +def test_element_summary_counts_and_concentration(): + per_species = { + # two species split at the SAME boundary with near-identical gaps + "A": [chunk(start=0, size=10, bi=1), chunk(start=43, size=5, bi=2)], + "B": [chunk(start=0, size=10, bi=1), chunk(start=44, size=5, bi=2)], + "C": [chunk(start=0, size=10, bi=1), chunk(start=10, size=5, bi=2)], + "D": [chunk(size=0, bi=1)], + "E": [chunk(start=0, size=10, bi=1)], + } + vals = dict(zip(ELEMENT_LOCI_HEADERS, summarizeElementLoci(per_species))) + assert vals["n.species"] == 5 + assert vals["n.split"] == 2 + assert vals["n.contiguous"] == 1 + assert vals["n.no.bases"] == 1 + assert vals["n.single"] == 1 + assert vals["split.bases.max"] == 34 # B's total + assert vals["split.gap.median"] == 33 # median of [33, 34] + assert vals["split.boundaries"] == 1 + assert vals["split.max.boundary.n.species"] == 2 + assert vals["split.max.boundary.gap.spread"] == 1 + + +def test_element_summary_separates_shared_from_diffuse(): + """ + max.boundary.n.species vs n.split is the shared-vs-independent signal: + equal means one event, much lower means scattered. + """ + shared = {f"sp{i}": [chunk(start=0, size=10, bi=1), chunk(start=43, size=5, bi=2)] for i in range(4)} + v = dict(zip(ELEMENT_LOCI_HEADERS, summarizeElementLoci(shared))) + assert v["n.split"] == 4 and v["split.max.boundary.n.species"] == 4 and v["split.boundaries"] == 1 + + diffuse = { + "a": [chunk(start=0, size=10, bi=1), chunk(start=43, size=5, bi=2)], + "b": [chunk(start=0, size=10, bi=2), chunk(start=43, size=5, bi=3)], + "c": [chunk(start=0, size=10, bi=3), chunk(start=43, size=5, bi=4)], + } + v = dict(zip(ELEMENT_LOCI_HEADERS, summarizeElementLoci(diffuse))) + assert v["n.split"] == 3 and v["split.max.boundary.n.species"] == 1 and v["split.boundaries"] == 3 + + +def test_summary_has_the_element_columns_and_matches_the_loci_table(tmp_path): + out = os.path.join(str(tmp_path), "out") + result = run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--loci-table"]) + assert result.returncode == 0, result.stderr + + summary = list(csv.DictReader(open(os.path.join(out, "maf_fetch_summary.tsv")), delimiter="\t")) + for col in ELEMENT_LOCI_HEADERS: + assert col in summary[0], col + + # the per-element counts must equal the loci table's own grouping + rows = read_loci(os.path.join(out, "maf_fetch_loci.tsv")) + seen, counts = set(), defaultdict(lambda: defaultdict(int)) + for r in rows: + key = (r["region_scaffold"], r["region_start"], r["region_end"]) + if (key, r["species"]) in seen: + continue + seen.add((key, r["species"])) + counts[key][r["class"]] += 1 + + checked = 0 + for r in summary: + key = (r["scaffold"], r["start"], r["end"]) + if key not in counts: + continue + checked += 1 + for cls, col in (("single", "n.single"), ("contiguous", "n.contiguous"), + ("split", "n.split"), ("multi_scaffold", "n.multi.scaffold"), + ("multi_strand", "n.multi.strand"), ("no_bases", "n.no.bases")): + assert counts[key][cls] == int(r[col]), (key, col) + assert checked > 0 + + +def test_element_columns_identical_in_maf_and_fasta_mode(tmp_path): + """ + The summary columns come from trimMafBlock's chunks on the MAF path and + from mafBlockToFasta's on the FASTA path; for a MAF with no duplicate + copies of a species per block they must agree exactly. + """ + maf_out = os.path.join(str(tmp_path), "maf") + fa_out = os.path.join(str(tmp_path), "fa") + common = [MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id"] + assert run_fetch(common + ["-o", maf_out]).returncode == 0 + assert run_fetch(common + ["-o", fa_out, "-f", "-fh", "species-coords-id"]).returncode == 0 + + def element_cols(path): + out = {} + for r in csv.DictReader(open(os.path.join(path, "maf_fetch_summary.tsv")), delimiter="\t"): + out[r["basename"]] = {c: r[c] for c in ELEMENT_LOCI_HEADERS} + return out + + assert element_cols(maf_out) == element_cols(fa_out) + + +@pytest.mark.parametrize("processes", [1, 4]) +def test_element_columns_deterministic_across_process_counts(processes, tmp_path): + base = os.path.join(str(tmp_path), "p1") + assert run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-o", base, "-p", "1"]).returncode == 0 + out = os.path.join(str(tmp_path), f"p{processes}") + assert run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-o", out, "-p", str(processes)]).returncode == 0 + + def rows(path): + return sorted( + tuple(r[c] for c in ["basename"] + ELEMENT_LOCI_HEADERS) + for r in csv.DictReader(open(os.path.join(path, "maf_fetch_summary.tsv")), delimiter="\t") + ) + assert rows(out) == rows(base) + + +# --------------------------------------------------------------------------- +# block_index / gap_to_next +# --------------------------------------------------------------------------- + + +def test_loci_table_block_index_and_gap_to_next(tmp_path): + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [("chr1", 24, 33, "cross")]) + out = os.path.join(str(tmp_path), "out") + assert run_fetch([MAF_FILE, bed, "--index", INDEX_FILE, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--loci-table"]).returncode == 0 + + rows = read_loci(os.path.join(out, "maf_fetch_loci.tsv")) + by_sp = defaultdict(list) + for r in rows: + by_sp[r["species"]].append(r) + + # human is split by 6bp between the region's two blocks; chimp is contiguous + human = sorted(by_sp["human"], key=lambda r: int(r["chunk_start"])) + assert [r["block_index"] for r in human] == ["1", "2"] + assert human[0]["gap_to_next"] == "6" # the gap sits at boundary 1->2 + assert human[1]["gap_to_next"] == "." # last chunk has no next + + chimp = sorted(by_sp["chimp"], key=lambda r: int(r["chunk_start"])) + assert chimp[0]["gap_to_next"] == "0" + + +def test_real_fixture_concentration_columns(tmp_path): + """ + tests/loci-classes.maf, region 3602239-3602958: 13 species split across 3 + boundaries with 13 sharing the busiest one -- a mostly-shared interruption. + """ + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [CLASSES_REGIONS[1]]) + out = os.path.join(str(tmp_path), "out") + assert run_fetch([CLASSES_MAF, bed, "--index", CLASSES_INDEX, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--fasta-dedupe", "most-seq"]).returncode == 0 + + r = next(csv.DictReader(open(os.path.join(out, "maf_fetch_summary.tsv")), delimiter="\t")) + assert int(r["n.split"]) == 13 + assert int(r["split.boundaries"]) == 3 + assert int(r["split.max.boundary.n.species"]) == 13 + assert int(r["split.max.boundary.gap.spread"]) == 1654 + assert int(r["n.multi.scaffold"]) == 1 + + +def test_no_bases_end_to_end_on_real_data(tmp_path): + """ + End-to-end cover for `no_bases`, the class most likely to silently corrupt + a downstream BED: filtering `single|contiguous` on the old behaviour would + have admitted these as clean, with zero-width intervals. + + A 5bp window at the very start of the fixture's first block leaves four + species with an s-line whose aligned bases all fall outside the trimmed + columns -- a real all-gap row, not a constructed one. + """ + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [("CM000994.3", 3541841, 3541846, "nobases")]) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([CLASSES_MAF, bed, "--index", CLASSES_INDEX, "-b", "id", "-o", out, + "-f", "-fh", "species-coords-id", "--fasta-dedupe", "most-seq", + "--loci-table"]) + assert result.returncode == 0, result.stderr + + rows = read_loci(os.path.join(out, "maf_fetch_loci.tsv")) + no_bases = {r["species"] for r in rows if r["class"] == "no_bases"} + assert len(no_bases) == 4, sorted(no_bases) + + # every no_bases row really is a zero-size chunk, and none is called clean + for r in rows: + if r["species"] in no_bases: + assert int(r["chunk_size"]) == 0, r + assert r["class"] not in ("single", "contiguous"), r + + # the summary count agrees with the loci table + summary = next(csv.DictReader(open(os.path.join(out, "maf_fetch_summary.tsv")), delimiter="\t")) + assert int(summary["n.no.bases"]) == 4 + assert int(summary["n.single"]) == 10 + + # and the FASTA rows are genuinely all-gap with a zero-width header span + seqs = {} + name = None + for line in open(os.path.join(out, "nobases.fa")): + if line.startswith(">"): + name = line.strip() + seqs[name] = "" + else: + seqs[name] += line.strip() + for header, seq in seqs.items(): + species = header[1:].split(".")[0] + if species in no_bases: + assert set(seq) == {"-"}, header + start, end = header.split(":")[1].split("(")[0].split("-") + assert start == end, header + + +# --------------------------------------------------------------------------- +# Reference-side summary columns +# --------------------------------------------------------------------------- + + +def test_summary_schema_marks_reference_columns(tmp_path): + out = os.path.join(str(tmp_path), "out") + assert run_fetch([MAF_FILE, BED_FILE, "--index", INDEX_FILE, "-b", "id", "-o", out]).returncode == 0 + + header = open(os.path.join(out, "maf_fetch_summary.tsv")).readline().rstrip("\n").split("\t") + assert "ref.n.overlapping.blocks" in header + assert "ref.block.bases" in header + # removed: it was structurally all-zeros and invited being read as a + # per-species contiguity check, which it never was + assert "interblock.distances" not in header + assert "block.lengths" not in header + + +def test_ref_coverage_gap_warns_and_shortens_ref_block_bases(tmp_path): + """ + 03-crossblocks-missing spans two blocks 6bp apart in the reference, so the + region holds reference positions no block covers. That used to surface only + as an `interblock.distances` value of [6]; it is now a warning, and + ref.block.bases falls short of the region width by exactly that much. + """ + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [("chr1", 24, 33, "cross")]) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([MAF_FILE, bed, "--index", INDEX_FILE, "-b", "id", "-o", out, "--verbose"]) + assert result.returncode == 0, result.stderr + combined = result.stdout + result.stderr + assert "ref-coverage-gap" in combined + assert "do not tile the reference" in combined + + r = next(csv.DictReader(open(os.path.join(out, "maf_fetch_summary.tsv")), delimiter="\t")) + width = int(r["end"]) - int(r["start"]) + assert int(r["ref.block.bases"]) == width - 6 + + +def test_no_ref_coverage_gap_warning_when_blocks_tile(tmp_path): + """ + A region inside a single block cannot have a reference coverage gap, so the + warning must not fire -- it would be noise on every normal run. + """ + bed = os.path.join(str(tmp_path), "r.bed") + write_bed(bed, [("chr2", 5, 10, "single")]) + out = os.path.join(str(tmp_path), "out") + + result = run_fetch([MAF_FILE, bed, "--index", INDEX_FILE, "-b", "id", "-o", out, "--verbose"]) + assert result.returncode == 0, result.stderr + assert "ref-coverage-gap" not in (result.stdout + result.stderr) + + r = next(csv.DictReader(open(os.path.join(out, "maf_fetch_summary.tsv")), delimiter="\t")) + assert int(r["ref.block.bases"]) == int(r["end"]) - int(r["start"])