"""Build a small, atomic SQLite taxonomy snapshot from NCBI metadata.""" import argparse from collections import Counter from datetime import datetime, timezone import hashlib import json from pathlib import Path import sqlite3 import tarfile import urllib.request ROOT = Path(__file__).resolve().parent SOURCES = { "assembly_summary_genbank.txt": "https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt", "taxdump.tar.gz": "https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz", } def dump_rows(archive, name): with archive.extractfile(name) as stream: for line in stream: yield [part.strip() for part in line.decode("utf-8").split("|")] def build(summary, taxdump, output, coverage_path=None): counts = Counter() covered_counts = Counter() inventory = json.loads(coverage_path.read_text()) if coverage_path and coverage_path.exists() else None covered = set(inventory["accessions"]) if inventory else set() matched = set() skipped = 0 with open(summary) as stream: for line in stream: if line.startswith("##"): continue if line.startswith("#"): header = line.lstrip("# ").rstrip("\n").split("\t") columns = {name: header.index(name) for name in ("assembly_accession", "taxid", "version_status")} continue row = line.rstrip("\n").split("\t") if not row[columns["assembly_accession"]].startswith("GCA_") or row[columns["version_status"]] != "latest": skipped += 1 continue try: taxid = int(row[columns["taxid"]]) except ValueError: taxid = -1 counts[taxid] += 1 accession = row[columns["assembly_accession"]] if accession in covered: if accession in matched: raise ValueError(f"Duplicate current assembly accession: {accession}") matched.add(accession) covered_counts[taxid] += 1 print(f"Read {sum(counts.values()):,} current assemblies", flush=True) parents, ranks, merged = {}, {}, {} with tarfile.open(taxdump, "r:gz") as archive: for row in dump_rows(archive, "nodes.dmp"): taxid = int(row[0]) parents[taxid], ranks[taxid] = int(row[1]), row[2] for row in dump_rows(archive, "merged.dmp"): merged[int(row[0])] = int(row[1]) direct = Counter() covered_direct = Counter() for taxid, count in counts.items(): original_taxid = taxid seen = set() while taxid in merged and taxid not in seen: seen.add(taxid) taxid = merged[taxid] direct[taxid if taxid in parents else -1] += count covered_direct[taxid if taxid in parents else -1] += covered_counts[original_taxid] parents[-1], ranks[-1] = 1, "unresolved" totals = Counter() covered_totals = Counter() for taxid, count in direct.items(): covered_count = covered_direct[taxid] seen = set() while True: if taxid in seen: raise ValueError("Taxonomy cycle detected") seen.add(taxid) totals[taxid] += count covered_totals[taxid] += covered_count if taxid == 1: break taxid = parents[taxid] names = {-1: "Unresolved taxonomy", 1: "GenBank genome assemblies"} for row in dump_rows(archive, "names.dmp"): taxid = int(row[0]) if taxid in totals and taxid != 1 and row[3] == "scientific name": names[taxid] = row[1] if not totals[1] or totals[1] != sum(counts.values()): raise ValueError("Invalid root assembly count") if covered_totals[1] != len(matched) or any(covered_totals[t] > totals[t] for t in totals): raise ValueError("Invalid covered assembly count") provenance = {} for path in (summary, taxdump): digest = hashlib.sha256() with open(path, "rb") as stream: for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): digest.update(chunk) source = {"url": SOURCES[path.name], "sha256": digest.hexdigest(), "bytes": path.stat().st_size} sidecar = path.with_suffix(path.suffix + ".source.json") if sidecar.exists(): source.update(json.loads(sidecar.read_text())) provenance[path.name] = source metadata = {"created_at": datetime.now(timezone.utc).isoformat(), "sources": provenance, "assembly_count": totals[1], "unresolved_assemblies": direct[-1], "taxa_with_assemblies": len(totals), "skipped_rows": skipped, "scope": "Current GCA assemblies in NCBI assembly_summary_genbank.txt; version_status=latest. Counts include all assembly levels and taxonomic groups. Not all GenBank nucleotide records."} if inventory: metadata["coverage"] = {k: v for k, v in inventory.items() if k not in ("accessions", "shards")} metadata["coverage"].update(assembly_versions_in_inventory=len(covered), matched_current_assemblies=len(matched), unmatched_assembly_versions=len(covered - matched), unmatched_accession_examples=sorted(covered - matched)[:25]) output.parent.mkdir(parents=True, exist_ok=True) temporary = output.with_suffix(".sqlite.tmp") temporary.unlink(missing_ok=True) try: with sqlite3.connect(temporary) as conn: conn.executescript("CREATE TABLE taxa (taxid INTEGER PRIMARY KEY, parent_id INTEGER, name TEXT, rank TEXT, direct_count INTEGER, total_count INTEGER, covered_direct INTEGER, covered_total INTEGER); CREATE INDEX parents ON taxa(parent_id, total_count DESC); CREATE TABLE metadata (value TEXT);") conn.executemany("INSERT INTO taxa VALUES (?,?,?,?,?,?,?,?)", ((t, parents[t], names.get(t, str(t)), ranks[t], direct[t], n, covered_direct[t], covered_totals[t]) for t, n in totals.items())) conn.execute("INSERT INTO metadata VALUES (?)", (json.dumps(metadata),)) temporary.replace(output) finally: temporary.unlink(missing_ok=True) print(json.dumps({k: v for k, v in metadata.items() if k != "sources"}, indent=2), flush=True) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--cache-dir", type=Path, default=ROOT / ".cache/coverage") parser.add_argument("--output", type=Path, default=ROOT / "data/taxonomy.sqlite") parser.add_argument("--download", action="store_true", help="Download fresh NCBI metadata before rebuilding (about 2 GB).") parser.add_argument("--coverage", type=Path, default=ROOT / "data/coverage.json", help="Published annotation inventory from refresh_coverage.py, used when present.") args = parser.parse_args() args.cache_dir.mkdir(parents=True, exist_ok=True) if args.download: for filename, url in SOURCES.items(): target = args.cache_dir / filename temporary = target.with_suffix(target.suffix + ".part") with urllib.request.urlopen(url, timeout=120) as response, open(temporary, "wb") as out: source = {"downloaded_at": datetime.now(timezone.utc).isoformat(), "last_modified": response.headers.get("Last-Modified")} while chunk := response.read(8 * 1024 * 1024): out.write(chunk) temporary.replace(target) target.with_suffix(target.suffix + ".source.json").write_text(json.dumps(source)) build(args.cache_dir / "assembly_summary_genbank.txt", args.cache_dir / "taxdump.tar.gz", args.output, args.coverage) if __name__ == "__main__": main()