File size: 7,962 Bytes
35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 7c1aeef 35d4af4 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 | """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()
|