From 578525bc951e1624887c6194c5d7b771c2132df9 Mon Sep 17 00:00:00 2001 From: Fang Chen Date: Thu, 16 Jul 2026 08:07:01 -0700 Subject: [PATCH] Phase 3: taxon milestone computation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add per-taxon assembly milestone tracking, computing when each taxon (at species and every higher rank) first reached four milestones — any assembly, EBP-affiliated assembly, meets-metric, and affiliated-and-meets-metric — from the version-tracked assembly TSVs that GoaT itself does not retain. - flows/lib/load_taxonomy.py: dev/test NCBI taxdump parser building the taxonomy contract (taxid -> scientific_name, rank, parent, canonical-rank lineage). Cycle-guarded parent-chain walks. Production uses upstream lineage. - flows/lib/compute_taxon_milestones.py: single chronological sweep producing species and higher-rank milestones; first-in-ranks recorded as a by-product; per-species latest_assembly / current / total counts. Emits compute.taxon_milestones.completed. Writes taxon_milestone_summary.tsv. - tests/test_taxon_milestones.py: 38 tests covering the parser, resolve_to_species, the sweep (each milestone earliest, tie-break, empty/None dates, in-ranks, latest-assembly), output schema, and validation invariants. - tests/README_phase_3_taxon_milestones.md: run/test guide and schema. --- flows/lib/compute_taxon_milestones.py | 412 +++++++++++++++++++ flows/lib/load_taxonomy.py | 180 ++++++++ tests/README_phase_3_taxon_milestones.md | 284 +++++++++++++ tests/test_taxon_milestones.py | 501 +++++++++++++++++++++++ 4 files changed, 1377 insertions(+) create mode 100644 flows/lib/compute_taxon_milestones.py create mode 100644 flows/lib/load_taxonomy.py create mode 100644 tests/README_phase_3_taxon_milestones.md create mode 100644 tests/test_taxon_milestones.py diff --git a/flows/lib/compute_taxon_milestones.py b/flows/lib/compute_taxon_milestones.py new file mode 100644 index 0000000..f34a473 --- /dev/null +++ b/flows/lib/compute_taxon_milestones.py @@ -0,0 +1,412 @@ +"""Compute per-taxon assembly milestones via a single chronological sweep. + +Reads assembly_current.tsv and assembly_historical.tsv (row-level version +detail) from the working directory, resolves each assembly's taxId to its +species, attaches the canonical-rank lineage, and walks all rows in +(releaseDate, accession) order to record, for every taxon at every rank, the +first assembly that reached each of four milestones: + + first_assembly — any assembly + first_ebp_assembly — EBP-affiliated submitter (PRJNA533106) + first_metric — meets the EBP quality standard (any submitter) + first_ebp_metric — meets the metric AND EBP-affiliated + +The any-submitter milestones (assembly, metric) additionally record, per +species, the canonical ranks at which that species was the first in its clade +to reach the milestone (first_assembly_in_ranks / first_metric_in_ranks). + +Writes taxon_milestone_summary.tsv (one row per taxon at any rank) and emits a +completion event so the rest of the pipeline can react. + +Taxonomy source: in production the lineage is attached to assembly rows +upstream; in dev/test this flow calls load_taxonomy.build_taxonomy on a local +NCBI taxdump (--taxdump_path). + +Usage: + SKIP_PREFECT=true python -m flows.lib.compute_taxon_milestones \\ + --work_dir tmp --taxdump_path test/taxonomy/ncbi +""" + +import csv +import os + +from flows.lib.conditional_import import emit_event, flow +from flows.lib.load_taxonomy import CANONICAL_RANKS, build_taxonomy +from flows.lib.shared_args import TAXDUMP_PATH, WORK_DIR +from flows.lib.shared_args import parse_args as _parse_args + +CURRENT_TSV = "assembly_current.tsv" +HISTORICAL_TSV = "assembly_historical.tsv" +OUTPUT_TSV = "taxon_milestone_summary.tsv" + +# EBP umbrella BioProject — affiliation is membership in this project. +EBP_BIOPROJECT = "PRJNA533106" + +# Levels walked per row, species first then up the canonical lineage. +LEVELS = ["species"] + ["genus", "family", "order", "class", "phylum", "kingdom"] + +# The four milestones. Each maps to its output column prefix; in_ranks is True +# for the any-submitter milestones that get a first_*_in_ranks column. +MILESTONES = [ + {"key": "assembly", "prefix": "first_assembly", "in_ranks": True}, + {"key": "ebp_assembly", "prefix": "first_ebp_assembly", "in_ranks": False}, + {"key": "metric", "prefix": "first_metric", "in_ranks": True}, + {"key": "ebp_metric", "prefix": "first_ebp_metric", "in_ranks": False}, +] + +OUTPUT_FIELDNAMES = [ + "taxid", + "rank", + "scientific_name", + "first_assembly_date", + "first_assembly_accession", + "first_ebp_assembly_date", + "first_ebp_assembly_accession", + "first_metric_date", + "first_metric_accession", + "first_ebp_metric_date", + "first_ebp_metric_accession", + "first_assembly_in_ranks", + "first_metric_in_ranks", + "latest_assembly_date", + "latest_assembly_accession", + "total_assemblies", + "current_assemblies", +] + + +def _accession(row: dict) -> str: + """Return the assembly accession, normalised across phases.""" + return row.get("accession") or row.get("genbankAccession") or "" + + +def _release_date(row: dict) -> str: + """Return the releaseDate, treating the literal 'None' as empty.""" + value = row.get("releaseDate") or "" + return "" if value == "None" else value + + +def _version_status(row: dict) -> str: + """Return versionStatus across camelCase / snake_case, default 'current'.""" + return row.get("versionStatus") or row.get("version_status") or "current" + + +def _is_affiliated(row: dict) -> bool: + """Return True if the row is under the EBP umbrella BioProject.""" + projects = (row.get("bioProjectAccession") or "").split(",") + return EBP_BIOPROJECT in (p.strip() for p in projects) + + +def _has_metric(row: dict) -> bool: + """Return True if the row meets the EBP quality standard.""" + value = row.get("ebpStandardDate") or "" + return value != "" and value != "None" + + +def milestone_predicate(row: dict, key: str) -> bool: + """Evaluate a milestone predicate against a row.""" + if key == "assembly": + return True + if key == "ebp_assembly": + return _is_affiliated(row) + if key == "metric": + return _has_metric(row) + if key == "ebp_metric": + return _has_metric(row) and _is_affiliated(row) + raise ValueError(f"unknown milestone key: {key}") + + +def load_assemblies(current_tsv: str, historical_tsv: str) -> list[dict]: + """Load and combine rows from the current and historical TSVs. + + Args: + current_tsv: Path to assembly_current.tsv. + historical_tsv: Path to assembly_historical.tsv. + + Returns: + Combined list of row dicts. + """ + rows = [] + for path, label in [(current_tsv, "current"), (historical_tsv, "historical")]: + if not os.path.exists(path): + print(f" Warning: {label} TSV not found: {path}") + continue + count = 0 + with open(path, encoding="utf-8") as f: + for row in csv.DictReader(f, delimiter="\t"): + rows.append(row) + count += 1 + print(f" Loaded {count} {label} rows") + return rows + + +def resolve_to_species(taxid: int, taxonomy: dict) -> int | None: + """Resolve a taxid to its species ancestor. + + Returns the taxid unchanged if it is species-rank, otherwise walks up the + parent chain to the first species ancestor. + + Args: + taxid: The taxid to resolve. + taxonomy: The taxonomy contract from build_taxonomy. + + Returns: + The species taxid, or None if no species ancestor exists. + """ + # A visited set guards against cycles in the parent chain (corruption or + # merged/deleted-taxid artifacts), which would otherwise loop forever. + visited = set() + node = taxonomy.get(taxid) + while node is not None and taxid not in visited: + if node["rank"] == "species": + return taxid + visited.add(taxid) + taxid = node.get("parent") + node = taxonomy.get(taxid) + return None + + +def compute_milestones(rows: list[dict], taxonomy: dict) -> tuple[dict, dict]: + """Run the single chronological sweep over all assembly rows. + + Args: + rows: Combined current + historical assembly rows. + taxonomy: The taxonomy contract. + + Returns: + A (taxa, species_extra) tuple: + taxa: {taxid: {rank, scientific_name, firsts: {milestone_key: + (date, accession)}}} for every taxon touched, at any rank. + species_extra: {species_taxid: {in_ranks: {milestone_key: [rank,...]}, + total_assemblies: int, current_assemblies: int, + latest: {date, accession} of the current/latest version}}. + """ + # Attach the resolved species taxid + lineage to each row; drop unresolvable. + enriched = [] + skipped = 0 + empty_date = 0 + for row in rows: + raw_taxid = row.get("taxId") or row.get("taxid") or "" + try: + taxid = int(raw_taxid) + except (ValueError, TypeError): + skipped += 1 + continue + species_taxid = resolve_to_species(taxid, taxonomy) + if species_taxid is None: + print(f" Skipping taxId {raw_taxid} ({_accession(row)}): no species ancestor") + skipped += 1 + continue + node = taxonomy[species_taxid] + enriched.append( + { + "row": row, + "species_taxid": species_taxid, + "lineage": node["lineage"], + "date": _release_date(row), + "accession": _accession(row), + } + ) + + # Per-species counts (include empty-date rows in totals) and the latest + # (current) assembly per species. + species_extra: dict[int, dict] = {} + for item in enriched: + sp = item["species_taxid"] + extra = species_extra.setdefault( + sp, + { + "in_ranks": {m["key"]: [] for m in MILESTONES if m["in_ranks"]}, + "total_assemblies": 0, + "current_assemblies": 0, + "latest": None, # (date, accession) of the current/latest version + }, + ) + extra["total_assemblies"] += 1 + is_current = _version_status(item["row"]) != "superseded" + if is_current: + extra["current_assemblies"] += 1 + # Latest = newest-dated current version; prefer current rows over + # superseded, then break ties by (date, accession). An empty date sorts + # first, so a dated current version always wins over an undated one. + candidate = (is_current, item["date"], item["accession"]) + if extra["latest"] is None or candidate > extra["latest"]["key"]: + extra["latest"] = { + "key": candidate, + "date": item["date"], + "accession": item["accession"], + } + + # Global chronological sort. Empty-date rows are excluded from the sweep + # (they cannot be a "first") but were already counted above. + dated = [item for item in enriched if item["date"]] + empty_date = len(enriched) - len(dated) + dated.sort(key=lambda item: (item["date"], item["accession"])) + + taxa: dict[int, dict] = {} + seen = {m["key"]: set() for m in MILESTONES} + + def ensure_taxon(taxid: int): + if taxid not in taxa: + node = taxonomy.get(taxid, {}) + taxa[taxid] = { + "rank": node.get("rank", ""), + "scientific_name": node.get("scientific_name", ""), + "firsts": {}, + } + + for item in dated: + sp = item["species_taxid"] + lineage = item["lineage"] + for milestone in MILESTONES: + key = milestone["key"] + if not milestone_predicate(item["row"], key): + continue + for level in LEVELS: + taxid = sp if level == "species" else lineage.get(level) + if not taxid or taxid in seen[key]: + continue + seen[key].add(taxid) + ensure_taxon(taxid) + taxa[taxid]["firsts"][key] = (item["date"], item["accession"]) + if milestone["in_ranks"] and level != "species": + species_extra[sp]["in_ranks"][key].append(level) + + # Every resolved species gets an output row even if all its rows had empty + # releaseDate (no milestone date, but its total/current counts still apply). + for sp in species_extra: + ensure_taxon(sp) + + print(f" Resolved {len(enriched)} rows ({skipped} skipped); {empty_date} with empty releaseDate excluded from dates") + print(f" Touched {len(taxa)} taxa across all ranks") + return taxa, species_extra + + +def build_output_rows(taxa: dict, species_extra: dict) -> list[dict]: + """Assemble the output rows from the sweep results. + + Args: + taxa: Per-taxon firsts from compute_milestones. + species_extra: Per-species in-ranks and counts. + + Returns: + List of output-row dicts keyed by OUTPUT_FIELDNAMES. + """ + out = [] + for taxid, info in taxa.items(): + is_species = info["rank"] == "species" + row = {name: "" for name in OUTPUT_FIELDNAMES} + row["taxid"] = taxid + row["rank"] = info["rank"] + row["scientific_name"] = info["scientific_name"] + + for milestone in MILESTONES: + first = info["firsts"].get(milestone["key"]) + date, accession = first if first else ("", "") + row[f"{milestone['prefix']}_date"] = date + row[f"{milestone['prefix']}_accession"] = accession + + if is_species: + extra = species_extra.get(taxid, {}) + in_ranks = extra.get("in_ranks", {}) + row["first_assembly_in_ranks"] = ",".join(in_ranks.get("assembly", [])) + row["first_metric_in_ranks"] = ",".join(in_ranks.get("metric", [])) + latest = extra.get("latest") + row["latest_assembly_date"] = latest["date"] if latest else "" + row["latest_assembly_accession"] = latest["accession"] if latest else "" + row["total_assemblies"] = extra.get("total_assemblies", 0) + row["current_assemblies"] = extra.get("current_assemblies", 0) + + out.append(row) + + # Stable, deterministic ordering: rank then taxid. + out.sort(key=lambda r: (str(r["rank"]), int(r["taxid"]))) + return out + + +@flow(log_prints=True) +def compute_taxon_milestones(work_dir: str = ".", taxdump_path: str | None = None) -> None: + """Compute per-taxon assembly milestones and write the summary TSV. + + Args: + work_dir: Directory containing assembly_current.tsv and + assembly_historical.tsv; output is written there too. + taxdump_path: Dev/test NCBI taxdump directory. In production the + lineage arrives via the upstream contract and this is omitted. + """ + current_tsv = os.path.join(work_dir, CURRENT_TSV) + historical_tsv = os.path.join(work_dir, HISTORICAL_TSV) + output_tsv = os.path.join(work_dir, OUTPUT_TSV) + + separator = "=" * 80 + print(f"\n{separator}") + print("TAXON MILESTONE COMPUTATION") + print(f"{separator}\n") + + print("[1/4] Loading taxonomy...") + if not taxdump_path: + # Production path would attach lineage upstream; without it there is + # nothing to resolve against, so fail loudly rather than silently. + raise SystemExit( + "No taxonomy source: pass --taxdump_path for the dev/test path, " + "or supply upstream lineage in production." + ) + taxonomy = build_taxonomy(taxdump_path) + + print("\n[2/4] Loading assemblies...") + rows = load_assemblies(current_tsv, historical_tsv) + + if not rows: + print(" No assembly data found — nothing to compute.") + emit_event( + event="compute.taxon_milestones.completed", + resource={ + "prefect.resource.id": f"compute.taxon_milestones.{work_dir}", + "prefect.resource.type": "taxon.milestones", + }, + payload={"taxa": 0, "status": "no_op"}, + ) + return + + print("\n[3/4] Computing milestones (single chronological sweep)...") + taxa, species_extra = compute_milestones(rows, taxonomy) + output_rows = build_output_rows(taxa, species_extra) + + print("\n[4/4] Writing output...") + with open(output_tsv, "w", newline="", encoding="utf-8") as f: + writer = csv.DictWriter(f, fieldnames=OUTPUT_FIELDNAMES, delimiter="\t") + writer.writeheader() + writer.writerows(output_rows) + + species_count = sum(1 for r in output_rows if r["rank"] == "species") + higher_count = len(output_rows) - species_count + + print(f"\n{separator}") + print("TAXON MILESTONE COMPUTATION COMPLETE") + print(f"{separator}") + print(f" Species rows: {species_count}") + print(f" Higher-rank rows: {higher_count}") + print(f" Output: {output_tsv}") + print(f"{separator}\n") + + emit_event( + event="compute.taxon_milestones.completed", + resource={ + "prefect.resource.id": f"compute.taxon_milestones.{work_dir}", + "prefect.resource.type": "taxon.milestones", + }, + payload={ + "taxa": len(output_rows), + "species": species_count, + "higher_rank": higher_count, + "status": "success", + }, + ) + + +if __name__ == "__main__": + args = _parse_args( + [WORK_DIR, TAXDUMP_PATH], + description="Compute per-taxon assembly milestones from version-tracked TSVs", + ) + compute_taxon_milestones(work_dir=args.work_dir, taxdump_path=args.taxdump_path) diff --git a/flows/lib/load_taxonomy.py b/flows/lib/load_taxonomy.py new file mode 100644 index 0000000..c8a4b42 --- /dev/null +++ b/flows/lib/load_taxonomy.py @@ -0,0 +1,180 @@ +"""Parse an NCBI taxdump into the Phase 3 taxonomy contract (dev/test only). + +This is a pure parser/builder. It performs no download and no caching — +fetching the taxdump is owned by update_ncbi_taxonomy.py. In production the +lineage is attached to assembly rows upstream (NCBI-dataset integration), so +this module is used only to run the pipeline end-to-end in development and in +tests. + +It reads the standard NCBI dump files: + + nodes.dmp: taxid | parent_taxid | rank | ... + names.dmp: taxid | name | unique_name | name_class | + +and builds the in-memory contract consumed by compute_taxon_milestones: + + { + taxid: { + 'scientific_name': str, + 'rank': str, # 'species', 'genus', 'subspecies', ... + 'parent': int, # one level up + 'lineage': {rank: taxid, ...} # canonical ranks only, genus..kingdom + }, + ... + } + +Usage: + python -m flows.lib.load_taxonomy --taxdump_path test/taxonomy/ncbi +""" + +import os + +from flows.lib.shared_args import TAXDUMP_PATH +from flows.lib.shared_args import parse_args as _parse_args + +# Canonical ranks captured in each node's lineage (genus up to kingdom). +CANONICAL_RANKS = {"genus", "family", "order", "class", "phylum", "kingdom"} + +NODES_FILE = "nodes.dmp" +NAMES_FILE = "names.dmp" + +# NCBI .dmp files separate fields with "\t|\t" and terminate rows with "\t|". +_FIELD_SEP = "\t|\t" +_ROW_TERMINATOR = "\t|" + + +def _split_dmp_line(line: str) -> list[str]: + """Split a single .dmp line into its fields. + + Strips the trailing row terminator ("\\t|") and the line break, then splits + on the field separator. Whitespace around each field is stripped. + + Args: + line: A raw line from a .dmp file. + + Returns: + List of field values (whitespace-trimmed). + """ + line = line.rstrip("\n").rstrip("\r") + if line.endswith(_ROW_TERMINATOR): + line = line[: -len(_ROW_TERMINATOR)] + return [field.strip() for field in line.split(_FIELD_SEP)] + + +def parse_nodes(nodes_dmp: str) -> dict[int, dict]: + """Parse nodes.dmp into a taxid -> {parent, rank} map. + + Args: + nodes_dmp: Path to the nodes.dmp file. + + Returns: + Dict mapping taxid (int) to {'parent': int, 'rank': str}. + """ + nodes = {} + with open(nodes_dmp, encoding="utf-8") as f: + for line in f: + if not line.strip(): + continue + fields = _split_dmp_line(line) + if len(fields) < 3: + continue + taxid = int(fields[0]) + parent = int(fields[1]) + rank = fields[2] + nodes[taxid] = {"parent": parent, "rank": rank} + return nodes + + +def parse_names(names_dmp: str) -> dict[int, str]: + """Parse names.dmp into a taxid -> scientific_name map. + + Only rows with name_class == 'scientific name' are retained. + + Args: + names_dmp: Path to the names.dmp file. + + Returns: + Dict mapping taxid (int) to its scientific name (str). + """ + names = {} + with open(names_dmp, encoding="utf-8") as f: + for line in f: + if not line.strip(): + continue + fields = _split_dmp_line(line) + if len(fields) < 4: + continue + if fields[3] == "scientific name": + names[int(fields[0])] = fields[1] + return names + + +def build_lineage(taxid: int, nodes: dict[int, dict]) -> dict[str, int]: + """Walk the parent chain and collect the taxid at each canonical rank. + + Starts from the node's parent and walks up to (but not including) the root + (taxid 1), recording the taxid wherever the node's rank is canonical. + + Args: + taxid: The taxid whose lineage to build. + nodes: The taxid -> {parent, rank} map from parse_nodes. + + Returns: + Dict mapping canonical rank name to the ancestor taxid at that rank. + """ + lineage = {} + node = nodes.get(taxid) + if node is None: + return lineage + # A visited set guards against cycles in the parent chain (corruption or + # merged/deleted-taxid artifacts), which would otherwise loop forever since + # the walk only stops at the root (taxid 1). + visited = {taxid} + current = node["parent"] + while current and current != 1 and current not in visited: + visited.add(current) + ancestor = nodes.get(current) + if ancestor is None: + break + if ancestor["rank"] in CANONICAL_RANKS: + lineage[ancestor["rank"]] = current + current = ancestor["parent"] + return lineage + + +def build_taxonomy(taxdump_dir: str) -> dict[int, dict]: + """Build the taxonomy contract from an NCBI taxdump directory. + + Combines nodes.dmp and names.dmp and computes a canonical-rank lineage for + every node. + + Args: + taxdump_dir: Directory containing nodes.dmp and names.dmp. + + Returns: + The taxonomy contract: taxid -> {scientific_name, rank, parent, lineage}. + """ + nodes_path = os.path.join(taxdump_dir, NODES_FILE) + names_path = os.path.join(taxdump_dir, NAMES_FILE) + + nodes = parse_nodes(nodes_path) + names = parse_names(names_path) + + taxonomy = {} + for taxid, node in nodes.items(): + taxonomy[taxid] = { + "scientific_name": names.get(taxid, ""), + "rank": node["rank"], + "parent": node["parent"], + "lineage": build_lineage(taxid, nodes), + } + + print(f" Loaded {len(taxonomy)} taxonomic nodes from {taxdump_dir}") + return taxonomy + + +if __name__ == "__main__": + args = _parse_args([TAXDUMP_PATH], description=__doc__) + if not args.taxdump_path: + raise SystemExit("--taxdump_path is required") + build_taxonomy(args.taxdump_path) diff --git a/tests/README_phase_3_taxon_milestones.md b/tests/README_phase_3_taxon_milestones.md new file mode 100644 index 0000000..eac73dd --- /dev/null +++ b/tests/README_phase_3_taxon_milestones.md @@ -0,0 +1,284 @@ +# Phase 3: Taxon Milestone Computation + +This guide covers the Phase 3 implementation, which computes per-taxon assembly +**milestone firsts** from the version-tracked assembly data produced by Phases +0–2 and writes `taxon_milestone_summary.tsv`. + +Phase 3 is **two modules**: + +- **`flows/lib/load_taxonomy.py`** (*dev/test only*) — parses a local NCBI + taxdump (`nodes.dmp` + `names.dmp`) into the taxonomy contract. In production + the lineage is attached to assembly rows upstream (NCBI-dataset integration), + so this module is used only to run the pipeline end-to-end without that join. +- **`flows/lib/compute_taxon_milestones.py`** — a single chronological sweep + that computes the two milestone deliverables — **species milestones** and + **higher-rank milestones** — then writes `taxon_milestone_summary.tsv`. The + **first-in-ranks** flags are not a third computation: they are a derived + species-level annotation, emitted as a by-product of the higher-rank step + (see below). + +Phase 0 (one-time backfill) and Phase 1 (daily incremental) must have run at +least once so that `assembly_current.tsv` and `assembly_historical.tsv` exist +before running Phase 3. Phase 2 and Phase 3 are **independent** — both trigger +from `update.assembly_versions.completed` and run in parallel; neither consumes +the other's output. + +## Files added in Phase 3 + +| File | Purpose | +|---|---| +| `flows/lib/load_taxonomy.py` | Pure NCBI-taxdump parser (no download, no cache); builds the taxonomy contract for the dev/test path | +| `flows/lib/compute_taxon_milestones.py` | Single-sweep milestone computation → `taxon_milestone_summary.tsv`; emits `compute.taxon_milestones.completed` | +| `tests/test_taxon_milestones.py` | Unit + end-to-end tests for both modules | +| `test/taxonomy/ncbi/` | Test taxdump fixture (`nodes.dmp` + `names.dmp`) for the *Isopoda* clade | + +## How it fits in the daily pipeline + +``` +parse_assembly_versions (writes assembly_historical.tsv; emits no event) + → update_assembly_versions → update.assembly_versions.completed + ├──► generate_assembly_summary → generate.assembly_summary.completed [Phase 2] + └──► compute_taxon_milestones → compute.taxon_milestones.completed [Phase 3] + (production: lineage attached upstream / supplied as the contract; + dev/test: load_taxonomy parses test/taxonomy/ncbi) +``` + +The trigger event `update.assembly_versions.completed` is emitted by +`update_assembly_versions.py` (the missing-version updater), **not** by +`parse_assembly_versions.py` (which writes the historical TSV but emits no +event). The updater emits the event on **every** run — including a `no_op` +emit when there are no missing versions to fetch — so Phase 2 and Phase 3 fan +out from it daily, in parallel, regardless of whether a backfill occurred. + +## The taxonomy contract + +Both environments feed the same in-memory structure; downstream code never +branches on the source: + +```python +CANONICAL_RANKS = {'genus', 'family', 'order', 'class', 'phylum', 'kingdom'} + +{ + taxid: { + 'scientific_name': str, + 'rank': str, # 'species', 'genus', 'subspecies', ... + 'parent': int, # one level up — used by resolve_to_species() + 'lineage': {rank: taxid, ...} # canonical ranks only, genus..kingdom + }, + ... +} +``` + +`scientific_name` and `parent` are present for **higher-rank** taxids too, not +just species — higher-rank output rows carry the rank taxon's own name, and +`resolve_to_species` walks `parent` to attribute subspecies to their species. + +`load_taxonomy.build_lineage` walks each node's parent chain up to the root +(taxid 1), recording the taxid at each canonical rank and skipping the +non-canonical intermediate ranks NCBI inserts (subphylum, superclass, suborder, +clade, …). + +## The four milestones + +All four are predicates over the **same** assembly rows on the **same** +`releaseDate` timeline: + +| Milestone | Date column (+`_accession`) | Predicate on a row | +|---|---|---| +| Any assembly | `first_assembly_date` | always true | +| EBP-affiliated assembly | `first_ebp_assembly_date` | `PRJNA533106` in `bioProjectAccession.split(',')` | +| Meets EBP metric (any submitter) | `first_metric_date` | `ebpStandardDate` non-empty | +| Meets metric **and** affiliated | `first_ebp_metric_date` | `ebpStandardDate` non-empty **AND** `PRJNA533106`-affiliated | + +**Naming convention (locked):** the `ebp_` prefix means *EBP-affiliated +submitter* (under the `PRJNA533106` umbrella project); `metric` means *meets the +EBP quality standard*, regardless of who submitted. + +> `ebpStandardDate` is set equal to a version's own `releaseDate` when EBP +> criteria pass (`check_ebp_criteria` in `flows/lib/utils.py`), so all four +> milestones share the `releaseDate` timeline. `first_ebp_metric_date` differs +> from `first_metric_date` only by the affiliation filter. +> +> Note on real data: `ebpStandardDate` and `releaseDate` may appear as the +> literal string `"None"` (not just empty) — Phase 3 treats `"None"` as empty +> for both the metric predicate and the milestone date. + +## The single chronological sweep + +``` +1. Load current + historical rows (csv module). +2. Per row: resolve taxId → species taxid (resolve_to_species); skip+log if + unresolvable. Attach species taxid + its canonical lineage. +3. Sort ALL dated rows by (releaseDate, accession) ascending — one global sort. + Rows with empty releaseDate are excluded from milestone dates but still + counted in total_assemblies; the count is logged. +4. For each milestone M, keep seen[M] = taxids already credited. Walk rows in + date order; for each row satisfying M, for each level in + [species, genus, family, order, class, phylum, kingdom]: + t = species_taxid (level==species) else lineage.get(rank) + if t and t not in seen[M]: + record (date, accession) as the first for (t, M); seen[M].add(t) + if M in {assembly, metric} and level != species: + append rank to that species's in_ranks[M] list +``` + +First-touch-in-date-order *is* that taxon's milestone — no aggregation, no +inverted index, no name-string matching. + +Note that the higher-rank milestone and the first-in-ranks flag are produced by +the **same** first-touch (step 4's inner block): the moment a higher-rank taxid +is first credited, the species that triggered it appends that rank to its +`in_ranks` list. So first-in-ranks is a re-attribution of the higher-rank +milestone back onto the species — not a separate pass. Only the two +any-submitter milestones (`assembly`, `metric`) do this re-attribution (the +`M in {assembly, metric}` guard); the two affiliated milestones still produce +higher-rank rows but no in-ranks flags. The earlier design computed first-in-ranks +in a separate third step that matched species to ranks by scientific-name string; +the single sweep removes that step entirely. + +Every resolved species also gets an output row even if all its assemblies have +empty dates (so its `total_assemblies` count is never lost). + +## Output: `taxon_milestone_summary.tsv` + +One row per taxon at any rank. Species rows carry in-ranks flags and counts; +higher-rank rows carry the earliest milestone date/accession for the clade. + +| Column | Species rows | Higher-rank rows | +|---|---|---| +| `taxid`, `rank`, `scientific_name` | ✓ | ✓ | +| `first_assembly_date` / `_accession` | ✓ | ✓ (earliest in clade) | +| `first_ebp_assembly_date` / `_accession` | ✓ | ✓ | +| `first_metric_date` / `_accession` | ✓ | ✓ | +| `first_ebp_metric_date` / `_accession` | ✓ | ✓ | +| `first_assembly_in_ranks` (e.g. `genus,family`) | ✓ | blank | +| `first_metric_in_ranks` | ✓ | blank | +| `latest_assembly_date` / `_accession` | ✓ | blank | +| `total_assemblies` | ✓ | blank | +| `current_assemblies` | ✓ | blank | + +`first_*_in_ranks` list the canonical ranks at which this species was the first +in its clade to reach that milestone (a subset of `CANONICAL_RANKS`). Only the +any-submitter milestones (`assembly`, `metric`) get in-ranks columns — +"first in clade to be assembled" and "first in clade to reach EBP quality" are +the scientifically meaningful firsts. + +`latest_assembly_date` / `latest_assembly_accession` report the species' +**current** (latest) version: the newest-dated non-superseded row, falling back +to the newest-dated row if none are current. `current_assemblies` counts rows +whose `versionStatus` (or `version_status`) is not `superseded`. + +## Assembly columns used + +| Column | Purpose | +|---|---| +| `taxId` | Taxonomy linkage (resolved to species) | +| `releaseDate` | Milestone date + global sort key | +| `bioProjectAccession` | EBP affiliation: `PRJNA533106` in `split(',')` | +| `ebpStandardDate` | Metric achieved when non-empty (and not `"None"`) | +| `genbankAccession` | Milestone accession + sort tiebreak (no `accession` column exists in the source TSVs) | +| `versionStatus` | current vs superseded (for `current_assemblies`) | + +## Prerequisites + +```bash +conda activate genomehubs_data +``` + +No network access or NCBI `datasets` CLI required — all inputs are local TSVs +plus the dev/test taxdump. + +## Running + +```bash +SKIP_PREFECT=true python3 -m flows.lib.compute_taxon_milestones \ + --work_dir /tmp/assembly-versions \ + --taxdump_path test/taxonomy/ncbi +``` + +Expects in `--work_dir`: +- `assembly_current.tsv` — written by `parse_ncbi_assemblies` (Phase 1 daily pipeline) +- `assembly_historical.tsv` — written by Phase 0 backfill, updated by Phase 1 incremental + +Writes to `--work_dir`: +- `taxon_milestone_summary.tsv` + +`--taxdump_path` is the dev/test taxonomy directory (containing `nodes.dmp` and +`names.dmp`). In production the lineage arrives via the upstream contract and +this argument is omitted; without either source the flow fails loudly rather +than producing an empty result. + +`load_taxonomy` can also be run standalone to inspect the parsed contract: + +```bash +SKIP_PREFECT=true python3 -m flows.lib.load_taxonomy \ + --taxdump_path test/taxonomy/ncbi +``` + +## Minimal end-to-end test + +```bash +mkdir -p /tmp/assembly-versions + +# Two species in different families of order Isopoda; Ceratothoa (2017) predates +# Asellus (2019), so it wins the higher-rank firsts at order and above. +cat > /tmp/assembly-versions/assembly_current.tsv << 'EOF' +genbankAccession taxId releaseDate bioProjectAccession ebpStandardDate versionStatus +GCA_AA.1 92525 2019-01-01 PRJNA533106 2019-01-01 current +EOF + +cat > /tmp/assembly-versions/assembly_historical.tsv << 'EOF' +genbankAccession taxId releaseDate bioProjectAccession ebpStandardDate versionStatus +GCA_CS.1 2922061 2017-01-01 PRJNA1 superseded +EOF + +SKIP_PREFECT=true python3 -m flows.lib.compute_taxon_milestones \ + --work_dir /tmp/assembly-versions \ + --taxdump_path test/taxonomy/ncbi + +cat /tmp/assembly-versions/taxon_milestone_summary.tsv +``` + +Illustrative Output: +- Order *Isopoda* (29979): `first_assembly_date=2017-01-01`, `accession=GCA_CS.1` + (earliest across the clade). +- *Asellus aquaticus* (92525): `first_assembly_in_ranks=genus,family` only — it + is first in its own genus/family but not at order and above (Ceratothoa took + those); `first_ebp_metric_date=2019-01-01` (affiliated + metric). +- *Ceratothoa steindachneri* (2922061): `first_assembly_in_ranks` = all six + canonical ranks; `current_assemblies=0` (its only version is superseded). + +## Validation invariants + +The predicates nest, so per taxon the milestone dates satisfy: + +``` +first_assembly_date ≤ first_metric_date ≤ first_ebp_metric_date +first_assembly_date ≤ first_ebp_assembly_date ≤ first_ebp_metric_date +``` + +(`first_metric_date` and `first_ebp_assembly_date` are not ordered relative to +each other.) Every `first_*_in_ranks` value is a subset of `CANONICAL_RANKS`. +These are asserted in the test suite. + +## Running the test suite + +```bash +SKIP_PREFECT=true python3 -m pytest tests/test_taxon_milestones.py -v +``` + +Covers: +- **load_taxonomy** — `parse_nodes`, `parse_names`, `build_lineage` (canonical + ranks only, species not in its own lineage), `build_taxonomy` against + `test/taxonomy/ncbi/`. +- **resolve_to_species** — species passthrough; subspecies → parent species; + genus / unknown taxid → `None`. +- **The sweep** — each of the four milestones picks the earliest qualifying row; + `(releaseDate, accession)` tie-breaking; the `"None"`-string trap; + empty-`releaseDate` excluded from dates but counted (incl. an empty-date-only + species still appearing in the output); in-ranks correctness across + genus/family/order; `current_assemblies` excludes superseded. +- **Output + invariants** — schema-exact rows; higher-rank rows blank their + in-ranks/counts; the nested date-ordering and `in_ranks ⊆ CANONICAL_RANKS` + invariants; an end-to-end flow run asserting the header and clade-earliest + dates. diff --git a/tests/test_taxon_milestones.py b/tests/test_taxon_milestones.py new file mode 100644 index 0000000..0b3a6db --- /dev/null +++ b/tests/test_taxon_milestones.py @@ -0,0 +1,501 @@ +"""Tests for load_taxonomy.py and compute_taxon_milestones.py (Phase 3). + +Covers: +- load_taxonomy: parse_nodes, parse_names, build_lineage, build_taxonomy + against the test/taxonomy/ncbi fixture. +- resolve_to_species: species passthrough, subspecies -> parent, unresolvable. +- The single chronological sweep: each of the four milestones picks the + earliest qualifying row; (releaseDate, accession) tie-breaking; empty + releaseDate excluded from dates but counted in total_assemblies; in-ranks + correctness across genus/family. +- Validation invariants: nested milestone-date ordering; in-ranks subset of + CANONICAL_RANKS. +""" + +import csv +import os +import sys +from pathlib import Path + +import pytest + +sys.path.insert(0, str(Path(__file__).parent.parent)) + +os.environ["SKIP_PREFECT"] = "true" + +from flows.lib.load_taxonomy import ( # noqa: E402 + CANONICAL_RANKS, + build_lineage, + build_taxonomy, + parse_names, + parse_nodes, +) +from flows.lib.compute_taxon_milestones import ( # noqa: E402 + OUTPUT_FIELDNAMES, + build_output_rows, + compute_milestones, + compute_taxon_milestones, + resolve_to_species, +) + +FIXTURE_DIR = str(Path(__file__).parent.parent / "test" / "taxonomy" / "ncbi") + +# Known taxids in the fixture (genus..kingdom for the Asellus aquaticus clade). +ASELLUS_AQUATICUS = 92525 # species +ASELLUS_GENUS = 92524 +ASELLIDAE_FAMILY = 63227 +ISOPODA_ORDER = 29979 +MALACOSTRACA_CLASS = 6681 +ARTHROPODA_PHYLUM = 6656 +METAZOA_KINGDOM = 33208 + +CERATOTHOA_STEINDACHNERI = 2922061 # species, different family/genus, same order +CERATOTHOA_GENUS = 432123 +CYMOTHOIDAE_FAMILY = 142082 + +JAERA_ISCHIOSETOSA = 2067965 # species +JAERA_PRAEHIRSUTA = 2613951 # species, same genus as above + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def write_tsv(path: Path, rows: list[dict]) -> None: + """Write a list of dicts to a tab-separated file.""" + if not rows: + path.write_text("", encoding="utf-8") + return + fieldnames = list(rows[0].keys()) + with open(path, "w", newline="", encoding="utf-8") as f: + writer = csv.DictWriter(f, fieldnames=fieldnames, delimiter="\t") + writer.writeheader() + writer.writerows(rows) + + +def read_tsv(path: Path) -> list[dict]: + """Read a tab-separated file into a list of dicts.""" + with open(path, encoding="utf-8") as f: + return list(csv.DictReader(f, delimiter="\t")) + + +def make_row(accession, taxid, release_date, bioproject="PRJNA1", ebp_date="", version_status="current"): + """Build a minimal assembly row dict.""" + return { + "genbankAccession": accession, + "taxId": str(taxid), + "releaseDate": release_date, + "bioProjectAccession": bioproject, + "ebpStandardDate": ebp_date, + "versionStatus": version_status, + } + + +@pytest.fixture(scope="module") +def taxonomy(): + """The taxonomy contract built from the test fixture.""" + return build_taxonomy(FIXTURE_DIR) + + +# --------------------------------------------------------------------------- +# load_taxonomy +# --------------------------------------------------------------------------- + +class TestParseNodes: + def test_parses_taxid_parent_rank(self): + nodes = parse_nodes(os.path.join(FIXTURE_DIR, "nodes.dmp")) + assert nodes[ASELLUS_AQUATICUS]["rank"] == "species" + assert nodes[ASELLUS_AQUATICUS]["parent"] == ASELLUS_GENUS + assert nodes[ARTHROPODA_PHYLUM]["rank"] == "phylum" + + def test_all_values_are_ints(self): + nodes = parse_nodes(os.path.join(FIXTURE_DIR, "nodes.dmp")) + for taxid, node in nodes.items(): + assert isinstance(taxid, int) + assert isinstance(node["parent"], int) + + +class TestParseNames: + def test_keeps_only_scientific_names(self): + names = parse_names(os.path.join(FIXTURE_DIR, "names.dmp")) + assert names[ASELLUS_AQUATICUS] == "Asellus aquaticus" + assert names[ARTHROPODA_PHYLUM] == "Arthropoda" + # 2759 has synonyms/common names but its scientific name is Eukaryota. + assert names[2759] == "Eukaryota" + + def test_synonyms_excluded(self): + names = parse_names(os.path.join(FIXTURE_DIR, "names.dmp")) + # 'Oniscus aquaticus' is a synonym of 92525, must not win. + assert names[ASELLUS_AQUATICUS] != "Oniscus aquaticus" + + +class TestBuildLineage: + def test_collects_canonical_ranks_only(self): + nodes = parse_nodes(os.path.join(FIXTURE_DIR, "nodes.dmp")) + lineage = build_lineage(ASELLUS_AQUATICUS, nodes) + assert lineage == { + "genus": ASELLUS_GENUS, + "family": ASELLIDAE_FAMILY, + "order": ISOPODA_ORDER, + "class": MALACOSTRACA_CLASS, + "phylum": ARTHROPODA_PHYLUM, + "kingdom": METAZOA_KINGDOM, + } + # Non-canonical ranks (subphylum, superclass, suborder...) excluded. + assert set(lineage.keys()) <= CANONICAL_RANKS + + def test_species_not_in_own_lineage(self): + nodes = parse_nodes(os.path.join(FIXTURE_DIR, "nodes.dmp")) + lineage = build_lineage(ASELLUS_AQUATICUS, nodes) + assert ASELLUS_AQUATICUS not in lineage.values() + + def test_cyclic_parent_chain_terminates(self): + # A 2-node cycle (and a self-parent) must not loop forever. + nodes = { + 10: {"parent": 20, "rank": "genus"}, + 20: {"parent": 10, "rank": "family"}, + 30: {"parent": 30, "rank": "order"}, # self-parent + } + assert build_lineage(10, nodes) == {"family": 20} + assert build_lineage(30, nodes) == {} + + +class TestBuildTaxonomy: + def test_contract_shape(self, taxonomy): + node = taxonomy[ASELLUS_AQUATICUS] + assert set(node.keys()) == {"scientific_name", "rank", "parent", "lineage"} + assert node["scientific_name"] == "Asellus aquaticus" + assert node["rank"] == "species" + assert node["parent"] == ASELLUS_GENUS + + def test_higher_rank_taxa_have_names(self, taxonomy): + assert taxonomy[ARTHROPODA_PHYLUM]["scientific_name"] == "Arthropoda" + assert taxonomy[ARTHROPODA_PHYLUM]["rank"] == "phylum" + + +# --------------------------------------------------------------------------- +# resolve_to_species +# --------------------------------------------------------------------------- + +class TestResolveToSpecies: + def test_species_passthrough(self, taxonomy): + assert resolve_to_species(ASELLUS_AQUATICUS, taxonomy) == ASELLUS_AQUATICUS + + def test_subspecies_resolves_to_parent_species(self, taxonomy): + # Inject a synthetic subspecies under Asellus aquaticus. + tax = dict(taxonomy) + tax[999001] = { + "scientific_name": "Asellus aquaticus subsp. test", + "rank": "subspecies", + "parent": ASELLUS_AQUATICUS, + "lineage": {}, + } + assert resolve_to_species(999001, tax) == ASELLUS_AQUATICUS + + def test_unresolvable_returns_none(self, taxonomy): + # Genus has no species ancestor. + assert resolve_to_species(ASELLUS_GENUS, taxonomy) is None + + def test_unknown_taxid_returns_none(self, taxonomy): + assert resolve_to_species(123456789, taxonomy) is None + + def test_cyclic_parent_chain_terminates(self): + # A cycle with no species in it must return None, not hang. + tax = { + 10: {"rank": "genus", "parent": 20, "scientific_name": "", "lineage": {}}, + 20: {"rank": "family", "parent": 10, "scientific_name": "", "lineage": {}}, + } + assert resolve_to_species(10, tax) is None + + def test_self_parent_terminates(self): + tax = {10: {"rank": "genus", "parent": 10, "scientific_name": "", "lineage": {}}} + assert resolve_to_species(10, tax) is None + + +# --------------------------------------------------------------------------- +# The chronological sweep +# --------------------------------------------------------------------------- + +class TestSweep: + def test_first_assembly_is_earliest_row(self, taxonomy): + rows = [ + make_row("GCA_002.1", ASELLUS_AQUATICUS, "2020-01-01"), + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-05-05"), + ] + taxa, _ = compute_milestones(rows, taxonomy) + assert taxa[ASELLUS_AQUATICUS]["firsts"]["assembly"] == ("2018-05-05", "GCA_001.1") + + def test_tiebreak_by_accession(self, taxonomy): + rows = [ + make_row("GCA_009.1", ASELLUS_AQUATICUS, "2019-01-01"), + make_row("GCA_003.1", ASELLUS_AQUATICUS, "2019-01-01"), + ] + taxa, _ = compute_milestones(rows, taxonomy) + assert taxa[ASELLUS_AQUATICUS]["firsts"]["assembly"] == ("2019-01-01", "GCA_003.1") + + def test_ebp_affiliation_predicate(self, taxonomy): + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-01-01", bioproject="PRJNA111"), + make_row("GCA_002.1", ASELLUS_AQUATICUS, "2019-01-01", bioproject="PRJNA533106,PRJNA222"), + ] + taxa, _ = compute_milestones(rows, taxonomy) + # first_assembly is the earliest regardless of affiliation. + assert taxa[ASELLUS_AQUATICUS]["firsts"]["assembly"][1] == "GCA_001.1" + # first_ebp_assembly is the affiliated one. + assert taxa[ASELLUS_AQUATICUS]["firsts"]["ebp_assembly"][1] == "GCA_002.1" + + def test_metric_predicate_and_none_string(self, taxonomy): + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-01-01", ebp_date="None"), + make_row("GCA_002.1", ASELLUS_AQUATICUS, "2019-01-01", ebp_date="2019-01-01"), + ] + taxa, _ = compute_milestones(rows, taxonomy) + # 'None' string must not count as metric; first metric is GCA_002.1. + assert taxa[ASELLUS_AQUATICUS]["firsts"]["metric"][1] == "GCA_002.1" + assert "metric" not in {k for k in taxa[ASELLUS_AQUATICUS]["firsts"] if False} + + def test_ebp_metric_requires_both(self, taxonomy): + rows = [ + # metric but not affiliated + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-01-01", bioproject="PRJNA111", ebp_date="2018-01-01"), + # affiliated but no metric + make_row("GCA_002.1", ASELLUS_AQUATICUS, "2019-01-01", bioproject="PRJNA533106", ebp_date=""), + # both + make_row("GCA_003.1", ASELLUS_AQUATICUS, "2020-01-01", bioproject="PRJNA533106", ebp_date="2020-01-01"), + ] + taxa, _ = compute_milestones(rows, taxonomy) + assert taxa[ASELLUS_AQUATICUS]["firsts"]["metric"][1] == "GCA_001.1" + assert taxa[ASELLUS_AQUATICUS]["firsts"]["ebp_assembly"][1] == "GCA_002.1" + assert taxa[ASELLUS_AQUATICUS]["firsts"]["ebp_metric"][1] == "GCA_003.1" + + def test_empty_date_excluded_but_counted(self, taxonomy): + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, ""), # no date + make_row("GCA_002.1", ASELLUS_AQUATICUS, "2019-01-01"), + ] + taxa, species_extra = compute_milestones(rows, taxonomy) + # First assembly skips the empty-date row. + assert taxa[ASELLUS_AQUATICUS]["firsts"]["assembly"][1] == "GCA_002.1" + # But it is still counted. + assert species_extra[ASELLUS_AQUATICUS]["total_assemblies"] == 2 + + def test_empty_date_only_species_still_in_output(self, taxonomy): + # A species whose ONLY assembly has an empty date must still appear as + # an output row with counts but no milestone dates. + rows = [make_row("GCA_001.1", ASELLUS_AQUATICUS, "")] + taxa, species_extra = compute_milestones(rows, taxonomy) + assert ASELLUS_AQUATICUS in taxa + assert taxa[ASELLUS_AQUATICUS]["firsts"] == {} + assert species_extra[ASELLUS_AQUATICUS]["total_assemblies"] == 1 + out = build_output_rows(taxa, species_extra) + sp_row = next(r for r in out if r["taxid"] == ASELLUS_AQUATICUS) + assert sp_row["first_assembly_date"] == "" + assert sp_row["total_assemblies"] == 1 + + def test_none_string_date_treated_as_empty(self, taxonomy): + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "None"), + make_row("GCA_002.1", ASELLUS_AQUATICUS, "2019-01-01"), + ] + taxa, _ = compute_milestones(rows, taxonomy) + assert taxa[ASELLUS_AQUATICUS]["firsts"]["assembly"][1] == "GCA_002.1" + + def test_current_assemblies_excludes_superseded(self, taxonomy): + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-01-01", version_status="superseded"), + make_row("GCA_001.2", ASELLUS_AQUATICUS, "2019-01-01", version_status="current"), + ] + _, species_extra = compute_milestones(rows, taxonomy) + assert species_extra[ASELLUS_AQUATICUS]["total_assemblies"] == 2 + assert species_extra[ASELLUS_AQUATICUS]["current_assemblies"] == 1 + + def test_latest_is_current_version(self, taxonomy): + # Latest = the current (non-superseded) version, even if an out-of-order + # superseded row has a newer date it must not win over a current one. + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-01-01", version_status="superseded"), + make_row("GCA_001.2", ASELLUS_AQUATICUS, "2021-01-01", version_status="current"), + ] + _, species_extra = compute_milestones(rows, taxonomy) + latest = species_extra[ASELLUS_AQUATICUS]["latest"] + assert (latest["date"], latest["accession"]) == ("2021-01-01", "GCA_001.2") + + def test_latest_prefers_current_over_newer_superseded(self, taxonomy): + # Pathological ordering: the superseded row is newer-dated than current. + # The current version must still be reported as latest. + rows = [ + make_row("GCA_001.2", ASELLUS_AQUATICUS, "2020-01-01", version_status="current"), + make_row("GCA_001.3", ASELLUS_AQUATICUS, "2022-01-01", version_status="superseded"), + ] + _, species_extra = compute_milestones(rows, taxonomy) + latest = species_extra[ASELLUS_AQUATICUS]["latest"] + assert latest["accession"] == "GCA_001.2" + + def test_latest_falls_back_when_all_superseded(self, taxonomy): + # If nothing is current, latest is the newest-dated superseded row. + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-01-01", version_status="superseded"), + make_row("GCA_001.2", ASELLUS_AQUATICUS, "2019-01-01", version_status="superseded"), + ] + _, species_extra = compute_milestones(rows, taxonomy) + latest = species_extra[ASELLUS_AQUATICUS]["latest"] + assert latest["accession"] == "GCA_001.2" + + def test_higher_rank_first_is_clade_earliest(self, taxonomy): + # Two species in the same order (Isopoda) but different families. + rows = [ + make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01"), + make_row("GCA_CS.1", CERATOTHOA_STEINDACHNERI, "2017-01-01"), + ] + taxa, _ = compute_milestones(rows, taxonomy) + # Order Isopoda's first assembly is the earliest across both clades. + assert taxa[ISOPODA_ORDER]["firsts"]["assembly"] == ("2017-01-01", "GCA_CS.1") + + def test_in_ranks_records_clade_firsts(self, taxonomy): + # Ceratothoa (2017) is earlier than Asellus (2019). Asellus is first in + # its OWN genus/family (different clade), Ceratothoa is first at order + # and above. So Asellus's in_ranks = genus, family (its clade) but NOT + # order/class/phylum/kingdom (Ceratothoa took those). + rows = [ + make_row("GCA_CS.1", CERATOTHOA_STEINDACHNERI, "2017-01-01"), + make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01"), + ] + _, species_extra = compute_milestones(rows, taxonomy) + asellus_in = species_extra[ASELLUS_AQUATICUS]["in_ranks"]["assembly"] + assert set(asellus_in) == {"genus", "family"} + cerat_in = species_extra[CERATOTHOA_STEINDACHNERI]["in_ranks"]["assembly"] + assert set(cerat_in) == {"genus", "family", "order", "class", "phylum", "kingdom"} + + def test_in_ranks_subset_of_canonical(self, taxonomy): + rows = [make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01")] + _, species_extra = compute_milestones(rows, taxonomy) + for milestone in ("assembly", "metric"): + ranks = species_extra[ASELLUS_AQUATICUS]["in_ranks"][milestone] + assert set(ranks) <= CANONICAL_RANKS + + +# --------------------------------------------------------------------------- +# Output rows + validation invariants +# --------------------------------------------------------------------------- + +class TestOutputRows: + def test_species_row_has_counts_and_in_ranks(self, taxonomy): + rows = [make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01")] + taxa, species_extra = compute_milestones(rows, taxonomy) + out = build_output_rows(taxa, species_extra) + species_row = next(r for r in out if r["taxid"] == ASELLUS_AQUATICUS) + assert species_row["rank"] == "species" + assert species_row["total_assemblies"] == 1 + assert species_row["first_assembly_in_ranks"] != "" + + def test_species_row_has_latest_assembly(self, taxonomy): + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2018-01-01", version_status="superseded"), + make_row("GCA_001.2", ASELLUS_AQUATICUS, "2021-05-05", version_status="current"), + ] + taxa, species_extra = compute_milestones(rows, taxonomy) + out = build_output_rows(taxa, species_extra) + species_row = next(r for r in out if r["taxid"] == ASELLUS_AQUATICUS) + assert species_row["latest_assembly_date"] == "2021-05-05" + assert species_row["latest_assembly_accession"] == "GCA_001.2" + + def test_higher_rank_row_blank_counts_and_in_ranks(self, taxonomy): + rows = [make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01")] + taxa, species_extra = compute_milestones(rows, taxonomy) + out = build_output_rows(taxa, species_extra) + phylum_row = next(r for r in out if r["taxid"] == ARTHROPODA_PHYLUM) + assert phylum_row["rank"] == "phylum" + assert phylum_row["total_assemblies"] == "" + assert phylum_row["current_assemblies"] == "" + assert phylum_row["latest_assembly_date"] == "" + assert phylum_row["latest_assembly_accession"] == "" + assert phylum_row["first_assembly_in_ranks"] == "" + # But it carries the clade's first assembly date. + assert phylum_row["first_assembly_date"] == "2019-01-01" + + def test_output_fieldnames_complete(self, taxonomy): + rows = [make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01")] + taxa, species_extra = compute_milestones(rows, taxonomy) + out = build_output_rows(taxa, species_extra) + for row in out: + assert set(row.keys()) == set(OUTPUT_FIELDNAMES) + + +class TestValidationInvariants: + def _dates(self, row): + return { + "assembly": row["first_assembly_date"], + "ebp_assembly": row["first_ebp_assembly_date"], + "metric": row["first_metric_date"], + "ebp_metric": row["first_ebp_metric_date"], + } + + def test_nested_date_ordering(self, taxonomy): + # A single species reaching all four milestones at distinct dates. + rows = [ + make_row("GCA_001.1", ASELLUS_AQUATICUS, "2015-01-01", bioproject="PRJNA111"), + make_row("GCA_002.1", ASELLUS_AQUATICUS, "2016-01-01", bioproject="PRJNA533106"), + make_row("GCA_003.1", ASELLUS_AQUATICUS, "2017-01-01", bioproject="PRJNA111", ebp_date="2017-01-01"), + make_row("GCA_004.1", ASELLUS_AQUATICUS, "2018-01-01", bioproject="PRJNA533106", ebp_date="2018-01-01"), + ] + taxa, species_extra = compute_milestones(rows, taxonomy) + out = build_output_rows(taxa, species_extra) + for row in out: + d = self._dates(row) + if d["assembly"] and d["metric"] and d["ebp_metric"]: + assert d["assembly"] <= d["metric"] <= d["ebp_metric"] + if d["assembly"] and d["ebp_assembly"] and d["ebp_metric"]: + assert d["assembly"] <= d["ebp_assembly"] <= d["ebp_metric"] + + def test_in_ranks_subset_invariant_all_rows(self, taxonomy): + rows = [ + make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01"), + make_row("GCA_CS.1", CERATOTHOA_STEINDACHNERI, "2017-01-01"), + ] + taxa, species_extra = compute_milestones(rows, taxonomy) + out = build_output_rows(taxa, species_extra) + for row in out: + for col in ("first_assembly_in_ranks", "first_metric_in_ranks"): + if row[col]: + assert set(row[col].split(",")) <= CANONICAL_RANKS + + +# --------------------------------------------------------------------------- +# End-to-end flow +# --------------------------------------------------------------------------- + +class TestFlow: + def test_end_to_end_writes_summary(self, tmp_path): + current = [ + make_row("GCA_AA.1", ASELLUS_AQUATICUS, "2019-01-01", bioproject="PRJNA533106", ebp_date="2019-01-01"), + make_row("GCA_JI.1", JAERA_ISCHIOSETOSA, "2020-06-06"), + ] + historical = [ + make_row("GCA_CS.1", CERATOTHOA_STEINDACHNERI, "2017-01-01", version_status="superseded"), + ] + write_tsv(tmp_path / "assembly_current.tsv", current) + write_tsv(tmp_path / "assembly_historical.tsv", historical) + + compute_taxon_milestones(work_dir=str(tmp_path), taxdump_path=FIXTURE_DIR) + + out = read_tsv(tmp_path / "taxon_milestone_summary.tsv") + assert out, "summary should not be empty" + by_taxid = {int(r["taxid"]): r for r in out} + + # Species rows present. + assert ASELLUS_AQUATICUS in by_taxid + assert by_taxid[ASELLUS_AQUATICUS]["rank"] == "species" + + # Order Isopoda's first assembly = earliest across the clade (Ceratothoa 2017). + assert by_taxid[ISOPODA_ORDER]["first_assembly_date"] == "2017-01-01" + assert by_taxid[ISOPODA_ORDER]["first_assembly_accession"] == "GCA_CS.1" + + # Header matches schema exactly. + with open(tmp_path / "taxon_milestone_summary.tsv", encoding="utf-8") as f: + header = f.readline().rstrip("\n").split("\t") + assert header == OUTPUT_FIELDNAMES + + def test_no_data_is_noop(self, tmp_path): + # No TSVs at all — flow should not raise and write no output. + compute_taxon_milestones(work_dir=str(tmp_path), taxdump_path=FIXTURE_DIR) + assert not (tmp_path / "taxon_milestone_summary.tsv").exists()