cgeorgiaw HF Staff commited on
Commit
7c1aeef
·
verified ·
1 Parent(s): 35d4af4

Color GenBank taxonomy by published annotation coverage

Browse files
FEATURE_BACKLOG.md CHANGED
@@ -8,8 +8,8 @@ Last updated: 2026-09-30.
8
 
9
  | ID | Feature | Status |
10
  | --- | --- | --- |
11
- | F001 | GenBank inference coverage map | Planned; not implemented |
12
- | F002 | Taxonomy map and navigation | GenBank assembly taxonomy implemented; annotation links and coverage pending |
13
  | F003 | Comparison with RefSeq annotations | Planned; not implemented |
14
  | F004 | Download annotations by accession, individually or in bulk | Partial; individual indexed segment downloads exist |
15
  | F005 | Fast accession lookup and retrieval across the full dataset | Prototype implemented for three remote objects; full coverage pending |
@@ -17,6 +17,10 @@ Last updated: 2026-09-30.
17
 
18
  ## F001 — GenBank inference coverage map
19
 
 
 
 
 
20
  **Goal:** Show how much of GenBank we currently have annotations for, and how that coverage grows as more inference completes.
21
 
22
  Desired behavior:
@@ -36,7 +40,7 @@ Decisions and dependencies for implementation:
36
 
37
  ## F002 — Taxonomy map and navigation
38
 
39
- Current implementation: the first tab shows an NCBI taxonomy sunburst sized by current GenBank assembly counts. It supports branch expansion, child/parent navigation, and scientific-name/taxon-ID search. A separate SQLite snapshot stores aggregated counts and provenance; `python refresh_taxonomy.py --download` refreshes it. See [refresh instructions](README.md#refresh-the-taxonomy-snapshot). The chart includes all assembly groups; it does not show inference coverage. Linking taxa to annotation accessions and coverage coloring remain future work. F001 is deferred for now.
40
 
41
  **Goal:** Provide a taxonomy-based view of the available genomes and annotations.
42
 
@@ -50,7 +54,7 @@ Desired behavior:
50
  Decisions and dependencies for implementation:
51
 
52
  - Implemented: sunburst with bounded subtrees and exact “Other taxa” counts; NCBI taxonomy with source checksums and dates; merged IDs resolved and unresolved IDs explicitly counted.
53
- - Remaining: connect the taxonomy inventory to inference coverage and accession-level annotation results.
54
 
55
  ## F003 — Comparison with RefSeq annotations
56
 
 
8
 
9
  | ID | Feature | Status |
10
  | --- | --- | --- |
11
+ | F001 | GenBank inference coverage map | Assembly-presence colors implemented; base completeness and history pending |
12
+ | F002 | Taxonomy map and navigation | Taxonomy and coverage colors implemented; annotation links pending |
13
  | F003 | Comparison with RefSeq annotations | Planned; not implemented |
14
  | F004 | Download annotations by accession, individually or in bulk | Partial; individual indexed segment downloads exist |
15
  | F005 | Fast accession lookup and retrieval across the full dataset | Prototype implemented for three remote objects; full coverage pending |
 
17
 
18
  ## F001 — GenBank inference coverage map
19
 
20
+ Current implementation: the taxonomy sunburst colors each group by the fraction of current GenBank assembly versions with published annotations. Hover shows counts and percentages; partial assemblies count as present, and this is explicitly distinguished from base coverage. The inventory includes completed published shards plus metadata from unmarked/legacy outputs, deduplicates assembly versions, and reports versions outside the current NCBI snapshot. `refresh_coverage.py` caches metadata by content fingerprints; `refresh_taxonomy.py` rolls coverage into SQLite. Refreshes are manual. See [coverage refresh and evidence limitations](README.md#refresh-annotation-coverage).
21
+
22
+ Remaining: partial versus complete assembly/base coverage, automatic updates and history, and navigation to available annotations. The existing completion markers do not bind input content hashes; stronger publication manifests would make the inventory independently auditable.
23
+
24
  **Goal:** Show how much of GenBank we currently have annotations for, and how that coverage grows as more inference completes.
25
 
26
  Desired behavior:
 
40
 
41
  ## F002 — Taxonomy map and navigation
42
 
43
+ Current implementation: the first tab shows an NCBI taxonomy sunburst sized by current GenBank assembly counts and colored by published annotation presence. It supports branch expansion, child/parent navigation, and scientific-name/taxon-ID search. A separate SQLite snapshot stores aggregated counts and provenance; see [refresh instructions](README.md#refresh-annotation-coverage). The chart includes all assembly groups. Linking taxa to annotation accessions remains future work.
44
 
45
  **Goal:** Provide a taxonomy-based view of the available genomes and annotations.
46
 
 
54
  Decisions and dependencies for implementation:
55
 
56
  - Implemented: sunburst with bounded subtrees and exact “Other taxa” counts; NCBI taxonomy with source checksums and dates; merged IDs resolved and unresolved IDs explicitly counted.
57
+ - Remaining: connect the taxonomy inventory to accession-level annotation results and distinguish partial/complete base coverage.
58
 
59
  ## F003 — Comparison with RefSeq annotations
60
 
README.md CHANGED
@@ -17,9 +17,9 @@ Future features and open design decisions are tracked in the [feature backlog](F
17
 
18
  ## GenBank taxonomy
19
 
20
- The first tab shows an interactive sunburst of **GenBank genome assemblies**, using NCBI's current GenBank assembly summary and taxonomy. Sector area represents assembly count. The second tab provides annotation lookup over the existing small index.
21
 
22
- The taxonomy includes all taxonomic groups and assembly levels in the current summary, accepting `GCA_` accessions with `version_status=latest`. Each assembly contributes once to its taxon and each ancestor. Merged taxon IDs are resolved; unresolved IDs are counted in a separate root branch. These counts describe genome assemblies, not all GenBank nucleotide records or annotation coverage.
23
 
24
  Click sectors to expand the loaded chart. Use the child-group menu or scientific-name/taxon-ID search to load another subtree, and the parent/home buttons to return. Each chart loads three descendant levels and the eight largest children per node; smaller children are combined into an accurately counted “Other taxa” sector. Search can reach those taxa individually. Hover percentages refer to the selected group; the group heading also shows its share of the full snapshot.
25
 
@@ -31,7 +31,22 @@ python refresh_taxonomy.py --download
31
 
32
  This downloads roughly 2 GB of **metadata only** to `.cache/coverage/`, then atomically rebuilds `data/taxonomy.sqlite`. To rebuild from already downloaded metadata, omit `--download`. No sequence or probability data is downloaded. The builder uses the Python standard library. The SQLite file contains taxon relationships, names, counts, and source URLs/checksums/dates; it is independent of `data/catalog.sqlite`.
33
 
34
- Upload the new `data/taxonomy.sqlite` to the Space and restart/redeploy to activate it. The app reads the bundled snapshot and does not download NCBI metadata at startup. Refreshes are manual; the displayed date is the snapshot build date. Annotation coverage coloring and accession navigation from taxonomy remain in the backlog.
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
35
 
36
  ## Run locally
37
 
 
17
 
18
  ## GenBank taxonomy
19
 
20
+ The first tab shows an interactive sunburst of **GenBank genome assemblies**, using NCBI's current GenBank assembly summary and taxonomy. Sector area represents assembly count. Color represents the fraction of assemblies with published annotations: gray for none found in the inventory, increasingly dark teal for higher percentages. Hover shows the numerator, denominator, and percentage. The second tab provides annotation lookup over the existing small index, which is independent of the broader coverage inventory.
21
 
22
+ The taxonomy includes all taxonomic groups and assembly levels in the current summary, accepting `GCA_` accessions with `version_status=latest`. Each assembly contributes once to its taxon and each ancestor. Merged taxon IDs are resolved; unresolved IDs are counted in a separate root branch. These counts describe genome assemblies, not all GenBank nucleotide records. An assembly counts as having annotations when it occurs in the publication inventory; this includes partial assemblies and does not imply full base coverage. Exact accession versions are matched to the NCBI snapshot. Published versions outside that snapshot are reported separately and excluded from chart percentages.
23
 
24
  Click sectors to expand the loaded chart. Use the child-group menu or scientific-name/taxon-ID search to load another subtree, and the parent/home buttons to return. Each chart loads three descendant levels and the eight largest children per node; smaller children are combined into an accurately counted “Other taxa” sector. Search can reach those taxa individually. Hover percentages refer to the selected group; the group heading also shows its share of the full snapshot.
25
 
 
31
 
32
  This downloads roughly 2 GB of **metadata only** to `.cache/coverage/`, then atomically rebuilds `data/taxonomy.sqlite`. To rebuild from already downloaded metadata, omit `--download`. No sequence or probability data is downloaded. The builder uses the Python standard library. The SQLite file contains taxon relationships, names, counts, and source URLs/checksums/dates; it is independent of `data/catalog.sqlite`.
33
 
34
+ Upload the new `data/taxonomy.sqlite` to the Space and restart/redeploy to activate it. The app reads the bundled snapshot and does not download NCBI metadata at startup. Refreshes are manual; the displayed dates distinguish the taxonomy build from the annotation inventory. Accession navigation from taxonomy and base-level completeness remain in the backlog.
35
+
36
+ ### Refresh annotation coverage
37
+
38
+ ```bash
39
+ python refresh_coverage.py
40
+ python refresh_taxonomy.py
41
+ ```
42
+
43
+ Use a Hugging Face token with read access to both `HuggingFaceBio/genbank-annotations` and `HuggingFaceBio/GENERanno-packed-inference-input`. The first command lists current published objects, identifies assemblies, and writes `data/coverage.json`. It reads Parquet footers and assembly-accession columns only; sequences and probabilities stay in the buckets. By default 12 workers process shards; change this with `--workers`.
44
+
45
+ For completed shards, the inventory joins `_SUCCESS.json` to assembly IDs in the corresponding packed `sources.parquet`. This relies on the pipeline's strict unpack contract: the publication script runs with `--strict`, checks part counts, and writes its completion marker after publishing annotations. It also assumes the input shard namespace still identifies the input used by that marker; the older markers do not contain input content hashes. This is assembly-presence evidence, not an audit of every output base. Shards without completion markers, including older single-file outputs, are inventoried directly from annotation metadata.
46
+
47
+ Assembly IDs are deduplicated across contigs, segments, and shards. The inventory records source/marker hashes and a fingerprint of published objects. Metadata results cached under `.cache/coverage/shards/` are reused only when those fingerprints match. A failed shard aborts publication of the new inventory; rerunning resumes using cached successes. Bucket listing is a dated observation, not an atomic bucket transaction.
48
+
49
+ The second command rebuilds taxonomy counts using that inventory and the cached NCBI files. Use `--download` to also refresh NCBI metadata. Upload both `data/coverage.json` (inventory and provenance) and `data/taxonomy.sqlite` (runtime counts) to the Space and redeploy. The refresh does not expand the accession lookup index.
50
 
51
  ## Run locally
52
 
data/coverage.json ADDED
The diff for this file is too large to render. See raw diff
 
data/taxonomy.sqlite CHANGED
@@ -1,3 +1,3 @@
1
  version https://git-lfs.github.com/spec/v1
2
- oid sha256:69388d96ee05e254b8eb8530a62a71d3d9e4f89e397427c1b7f99c140d488b55
3
- size 24084480
 
1
  version https://git-lfs.github.com/spec/v1
2
+ oid sha256:3bc92845313c3bfc3f406a9ded21d32ebc3dcc336d44a1eb791ea2281094e62c
3
+ size 24797184
refresh_coverage.py ADDED
@@ -0,0 +1,175 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ """Inventory published annotation assemblies using completion markers and metadata.
2
+
3
+ Completed shards use assembly IDs from their packed sources, under the pipeline's
4
+ strict unpack/publication contract. Shards without markers read annotation IDs
5
+ directly. No sequences or probability columns are read. Cached results are tied
6
+ to source and publication content hashes, so reruns reuse unchanged shards.
7
+ """
8
+ import argparse
9
+ from collections import defaultdict
10
+ from concurrent.futures import ThreadPoolExecutor, as_completed
11
+ from datetime import datetime, timezone
12
+ import hashlib
13
+ import json
14
+ from pathlib import Path
15
+ import re
16
+ import time
17
+
18
+ from huggingface_hub import HfApi, HfFileSystem, get_token
19
+ import pyarrow.parquet as pq
20
+
21
+ from remote_catalog import MeasuredFile, source_info
22
+
23
+ ROOT = Path(__file__).resolve().parent
24
+ ANNOTATIONS = "HuggingFaceBio/genbank-annotations"
25
+ INPUTS = "HuggingFaceBio/GENERanno-packed-inference-input"
26
+
27
+
28
+ def fingerprint(value):
29
+ return hashlib.sha256(json.dumps(value, sort_keys=True).encode()).hexdigest()
30
+
31
+
32
+ def read_accessions(fs, bucket, path):
33
+ accessions = set()
34
+ with MeasuredFile(fs, f"buckets/{bucket}/{path}") as stream:
35
+ pf = pq.ParquetFile(stream)
36
+ column = pf.schema.names.index("assembly_accession")
37
+ for group in range(pf.num_row_groups):
38
+ stats = pf.metadata.row_group(group).column(column).statistics
39
+ if stats and stats.has_min_max and stats.min == stats.max and stats.null_count == 0:
40
+ accessions.add(stats.min)
41
+ else:
42
+ accessions.update(pf.read_row_group(group, columns=["assembly_accession"], use_threads=False)["assembly_accession"].to_pylist())
43
+ fetched = stream.bytes_read
44
+ if any(not isinstance(a, str) or not re.fullmatch(r"GCA_\d+\.\d+", a) for a in accessions):
45
+ raise ValueError(f"Invalid assembly accession in {path}")
46
+ return accessions, fetched
47
+
48
+
49
+ def collect(cache_dir, output, workers):
50
+ cache_dir.mkdir(parents=True, exist_ok=True)
51
+ token = get_token()
52
+ api = HfApi(token=token)
53
+ started = datetime.now(timezone.utc).isoformat()
54
+ shards = defaultdict(lambda: {"objects": [], "marker": None})
55
+ print("Listing published annotation objects…", flush=True)
56
+ for item in api.list_bucket_tree(ANNOTATIONS, "annotations", recursive=True):
57
+ if item.type != "file":
58
+ continue
59
+ parts = item.path.split("/")
60
+ if len(parts) == 3 and parts[-1].endswith(".parquet"):
61
+ # Older publishers wrote one file per shard, without a shard folder.
62
+ shards["/".join(parts[1:])]["objects"].append((item.path, item.xet_hash, item.size))
63
+ continue
64
+ if len(parts) != 4:
65
+ if item.path.endswith(".parquet"):
66
+ raise ValueError(f"Unrecognized annotation layout: {item.path}")
67
+ continue
68
+ shard = "/".join(parts[1:3])
69
+ if parts[-1].endswith(".parquet"):
70
+ shards[shard]["objects"].append((item.path, item.xet_hash, item.size))
71
+ elif parts[-1] == "_SUCCESS.json":
72
+ shards[shard]["marker"] = (item.path, item.xet_hash)
73
+ shards = {key: val for key, val in shards.items() if val["objects"]}
74
+ if not shards:
75
+ raise ValueError("No annotation objects found")
76
+ print(f"Found {len(shards):,} published shards, {sum(len(v['objects']) for v in shards.values()):,} objects", flush=True)
77
+
78
+ def scan(shard, entry):
79
+ fs, worker_api = HfFileSystem(token=token), HfApi(token=token)
80
+ source = None
81
+ if entry["marker"]:
82
+ path = f"final_overlapped_pack/packed_blocks_by_division/{shard}/sources.parquet"
83
+ info = source_info(worker_api, INPUTS, path)
84
+ if info is None:
85
+ raise ValueError(f"Missing packed source metadata for {shard}")
86
+ source = (path, info.xet_hash)
87
+ identity = {"objects": sorted(entry["objects"]), "marker": entry["marker"], "source": source, "version": 1}
88
+ digest = fingerprint(identity)
89
+ target = cache_dir / (shard.replace("/", "__") + ".json")
90
+ if target.exists():
91
+ previous = json.loads(target.read_text())
92
+ if previous["fingerprint"] == digest:
93
+ return previous
94
+ accessions, fetched = set(), 0
95
+ marker = None
96
+ if source:
97
+ marker = json.loads(fs.cat(f"buckets/{ANNOTATIONS}/{entry['marker'][0]}"))
98
+ division, shard_id = shard.split("/")
99
+ if (marker.get("division") != division or marker.get("shard") != shard_id
100
+ or int(marker.get("expected_parts", 0)) <= 0
101
+ or int(marker.get("pack_count", 0)) < int(marker["expected_parts"])
102
+ or int(marker.get("prob_count", 0)) < int(marker["expected_parts"])):
103
+ raise ValueError(f"Invalid completion marker for {shard}")
104
+ accessions, fetched = read_accessions(fs, INPUTS, source[0])
105
+ after = source_info(worker_api, INPUTS, source[0])
106
+ if after is None or after.xet_hash != source[1]:
107
+ raise ValueError(f"Source changed while reading {shard}; retry")
108
+ after_marker = source_info(worker_api, ANNOTATIONS, entry['marker'][0])
109
+ if after_marker is None or after_marker.xet_hash != entry['marker'][1]:
110
+ raise ValueError(f"Completion marker changed while reading {shard}; retry")
111
+ else:
112
+ for path, content_hash, _ in entry["objects"]:
113
+ ids, size = read_accessions(fs, ANNOTATIONS, path)
114
+ accessions.update(ids)
115
+ fetched += size
116
+ after = source_info(worker_api, ANNOTATIONS, path)
117
+ if after is None or after.xet_hash != content_hash:
118
+ raise ValueError(f"Annotation changed while reading {path}; retry")
119
+ result = {"shard": shard, "fingerprint": digest, "accessions": sorted(accessions),
120
+ "annotation_objects": len(entry["objects"]), "metadata_bytes_read": fetched,
121
+ "method": "strict_completion_marker_and_packed_sources" if source else "annotation_metadata",
122
+ "source": source, "marker": entry["marker"], "completion": marker}
123
+ temporary = target.with_suffix(".tmp")
124
+ temporary.write_text(json.dumps(result))
125
+ temporary.replace(target)
126
+ return result
127
+
128
+ def retry_scan(shard, entry):
129
+ for attempt in range(3):
130
+ try:
131
+ return scan(shard, entry)
132
+ except Exception:
133
+ if attempt == 2:
134
+ raise
135
+ time.sleep(2 ** attempt)
136
+
137
+ results, errors = [], []
138
+ with ThreadPoolExecutor(max_workers=workers) as pool:
139
+ futures = {pool.submit(retry_scan, shard, entry): shard for shard, entry in shards.items()}
140
+ for future in as_completed(futures):
141
+ try:
142
+ results.append(future.result())
143
+ except Exception as exc:
144
+ errors.append(f"{futures[future]}: {type(exc).__name__}: {exc}")
145
+ print(errors[-1], flush=True)
146
+ if (len(results) + len(errors)) % 20 == 0:
147
+ print(f"Inventoried {len(results):,}/{len(shards):,} shards; {len(errors)} errors", flush=True)
148
+ if errors:
149
+ raise RuntimeError(f"Coverage snapshot not replaced: {len(errors)} shards failed. Retry to reuse cached successes.\n" + "\n".join(errors))
150
+ accessions = sorted({a for r in results for a in r["accessions"]})
151
+ inventory = {"created_at": datetime.now(timezone.utc).isoformat(), "listing_started_at": started,
152
+ "annotation_bucket": ANNOTATIONS, "input_bucket": INPUTS,
153
+ "scope": "Assemblies with published annotations; includes partial assemblies. Exact accession versions; no base completeness claim.",
154
+ "method": "Completed shard assembly IDs from packed sources under strict unpack publication contract; unmarked shards read from annotation metadata. Source namespace is assumed to remain associated with its published completion marker.",
155
+ "annotation_object_count": sum(r["annotation_objects"] for r in results),
156
+ "shard_count": len(results), "metadata_bytes_read": sum(r["metadata_bytes_read"] for r in results),
157
+ "inventory_fingerprint": fingerprint(sorted(r["fingerprint"] for r in results)),
158
+ "accessions": accessions,
159
+ "shards": [{k: v for k, v in r.items() if k != "accessions"} for r in sorted(results, key=lambda r: r['shard'])]}
160
+ output.parent.mkdir(parents=True, exist_ok=True)
161
+ temporary = output.with_suffix(".json.tmp")
162
+ temporary.write_text(json.dumps(inventory))
163
+ temporary.replace(output)
164
+ print(f"Saved {len(accessions):,} distinct assembly versions to {output}", flush=True)
165
+
166
+
167
+ if __name__ == "__main__":
168
+ parser = argparse.ArgumentParser(description=__doc__)
169
+ parser.add_argument("--cache-dir", type=Path, default=ROOT / ".cache/coverage/shards")
170
+ parser.add_argument("--output", type=Path, default=ROOT / "data/coverage.json")
171
+ parser.add_argument("--workers", type=int, default=12)
172
+ args = parser.parse_args()
173
+ if not 1 <= args.workers <= 32:
174
+ parser.error("workers must be between 1 and 32")
175
+ collect(args.cache_dir, args.output, args.workers)
refresh_taxonomy.py CHANGED
@@ -22,8 +22,12 @@ def dump_rows(archive, name):
22
  yield [part.strip() for part in line.decode("utf-8").split("|")]
23
 
24
 
25
- def build(summary, taxdump, output):
26
  counts = Counter()
 
 
 
 
27
  skipped = 0
28
  with open(summary) as stream:
29
  for line in stream:
@@ -42,6 +46,12 @@ def build(summary, taxdump, output):
42
  except ValueError:
43
  taxid = -1
44
  counts[taxid] += 1
 
 
 
 
 
 
45
  print(f"Read {sum(counts.values()):,} current assemblies", flush=True)
46
  parents, ranks, merged = {}, {}, {}
47
  with tarfile.open(taxdump, "r:gz") as archive:
@@ -51,21 +61,27 @@ def build(summary, taxdump, output):
51
  for row in dump_rows(archive, "merged.dmp"):
52
  merged[int(row[0])] = int(row[1])
53
  direct = Counter()
 
54
  for taxid, count in counts.items():
 
55
  seen = set()
56
  while taxid in merged and taxid not in seen:
57
  seen.add(taxid)
58
  taxid = merged[taxid]
59
  direct[taxid if taxid in parents else -1] += count
 
60
  parents[-1], ranks[-1] = 1, "unresolved"
61
  totals = Counter()
 
62
  for taxid, count in direct.items():
 
63
  seen = set()
64
  while True:
65
  if taxid in seen:
66
  raise ValueError("Taxonomy cycle detected")
67
  seen.add(taxid)
68
  totals[taxid] += count
 
69
  if taxid == 1:
70
  break
71
  taxid = parents[taxid]
@@ -76,6 +92,8 @@ def build(summary, taxdump, output):
76
  names[taxid] = row[1]
77
  if not totals[1] or totals[1] != sum(counts.values()):
78
  raise ValueError("Invalid root assembly count")
 
 
79
  provenance = {}
80
  for path in (summary, taxdump):
81
  digest = hashlib.sha256()
@@ -91,18 +109,24 @@ def build(summary, taxdump, output):
91
  "assembly_count": totals[1], "unresolved_assemblies": direct[-1],
92
  "taxa_with_assemblies": len(totals), "skipped_rows": skipped,
93
  "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."}
 
 
 
 
 
 
94
  output.parent.mkdir(parents=True, exist_ok=True)
95
  temporary = output.with_suffix(".sqlite.tmp")
96
  temporary.unlink(missing_ok=True)
97
  try:
98
  with sqlite3.connect(temporary) as conn:
99
- conn.executescript("CREATE TABLE taxa (taxid INTEGER PRIMARY KEY, parent_id INTEGER, name TEXT, rank TEXT, direct_count INTEGER, total_count INTEGER); CREATE INDEX parents ON taxa(parent_id, total_count DESC); CREATE TABLE metadata (value TEXT);")
100
- conn.executemany("INSERT INTO taxa VALUES (?,?,?,?,?,?)", ((t, parents[t], names.get(t, str(t)), ranks[t], direct[t], n) for t, n in totals.items()))
101
  conn.execute("INSERT INTO metadata VALUES (?)", (json.dumps(metadata),))
102
  temporary.replace(output)
103
  finally:
104
  temporary.unlink(missing_ok=True)
105
- print(json.dumps(metadata, indent=2), flush=True)
106
 
107
 
108
  def main():
@@ -110,6 +134,7 @@ def main():
110
  parser.add_argument("--cache-dir", type=Path, default=ROOT / ".cache/coverage")
111
  parser.add_argument("--output", type=Path, default=ROOT / "data/taxonomy.sqlite")
112
  parser.add_argument("--download", action="store_true", help="Download fresh NCBI metadata before rebuilding (about 2 GB).")
 
113
  args = parser.parse_args()
114
  args.cache_dir.mkdir(parents=True, exist_ok=True)
115
  if args.download:
@@ -122,7 +147,7 @@ def main():
122
  out.write(chunk)
123
  temporary.replace(target)
124
  target.with_suffix(target.suffix + ".source.json").write_text(json.dumps(source))
125
- build(args.cache_dir / "assembly_summary_genbank.txt", args.cache_dir / "taxdump.tar.gz", args.output)
126
 
127
 
128
  if __name__ == "__main__":
 
22
  yield [part.strip() for part in line.decode("utf-8").split("|")]
23
 
24
 
25
+ def build(summary, taxdump, output, coverage_path=None):
26
  counts = Counter()
27
+ covered_counts = Counter()
28
+ inventory = json.loads(coverage_path.read_text()) if coverage_path and coverage_path.exists() else None
29
+ covered = set(inventory["accessions"]) if inventory else set()
30
+ matched = set()
31
  skipped = 0
32
  with open(summary) as stream:
33
  for line in stream:
 
46
  except ValueError:
47
  taxid = -1
48
  counts[taxid] += 1
49
+ accession = row[columns["assembly_accession"]]
50
+ if accession in covered:
51
+ if accession in matched:
52
+ raise ValueError(f"Duplicate current assembly accession: {accession}")
53
+ matched.add(accession)
54
+ covered_counts[taxid] += 1
55
  print(f"Read {sum(counts.values()):,} current assemblies", flush=True)
56
  parents, ranks, merged = {}, {}, {}
57
  with tarfile.open(taxdump, "r:gz") as archive:
 
61
  for row in dump_rows(archive, "merged.dmp"):
62
  merged[int(row[0])] = int(row[1])
63
  direct = Counter()
64
+ covered_direct = Counter()
65
  for taxid, count in counts.items():
66
+ original_taxid = taxid
67
  seen = set()
68
  while taxid in merged and taxid not in seen:
69
  seen.add(taxid)
70
  taxid = merged[taxid]
71
  direct[taxid if taxid in parents else -1] += count
72
+ covered_direct[taxid if taxid in parents else -1] += covered_counts[original_taxid]
73
  parents[-1], ranks[-1] = 1, "unresolved"
74
  totals = Counter()
75
+ covered_totals = Counter()
76
  for taxid, count in direct.items():
77
+ covered_count = covered_direct[taxid]
78
  seen = set()
79
  while True:
80
  if taxid in seen:
81
  raise ValueError("Taxonomy cycle detected")
82
  seen.add(taxid)
83
  totals[taxid] += count
84
+ covered_totals[taxid] += covered_count
85
  if taxid == 1:
86
  break
87
  taxid = parents[taxid]
 
92
  names[taxid] = row[1]
93
  if not totals[1] or totals[1] != sum(counts.values()):
94
  raise ValueError("Invalid root assembly count")
95
+ if covered_totals[1] != len(matched) or any(covered_totals[t] > totals[t] for t in totals):
96
+ raise ValueError("Invalid covered assembly count")
97
  provenance = {}
98
  for path in (summary, taxdump):
99
  digest = hashlib.sha256()
 
109
  "assembly_count": totals[1], "unresolved_assemblies": direct[-1],
110
  "taxa_with_assemblies": len(totals), "skipped_rows": skipped,
111
  "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."}
112
+ if inventory:
113
+ metadata["coverage"] = {k: v for k, v in inventory.items() if k not in ("accessions", "shards")}
114
+ metadata["coverage"].update(assembly_versions_in_inventory=len(covered),
115
+ matched_current_assemblies=len(matched),
116
+ unmatched_assembly_versions=len(covered - matched),
117
+ unmatched_accession_examples=sorted(covered - matched)[:25])
118
  output.parent.mkdir(parents=True, exist_ok=True)
119
  temporary = output.with_suffix(".sqlite.tmp")
120
  temporary.unlink(missing_ok=True)
121
  try:
122
  with sqlite3.connect(temporary) as conn:
123
+ 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);")
124
+ 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()))
125
  conn.execute("INSERT INTO metadata VALUES (?)", (json.dumps(metadata),))
126
  temporary.replace(output)
127
  finally:
128
  temporary.unlink(missing_ok=True)
129
+ print(json.dumps({k: v for k, v in metadata.items() if k != "sources"}, indent=2), flush=True)
130
 
131
 
132
  def main():
 
134
  parser.add_argument("--cache-dir", type=Path, default=ROOT / ".cache/coverage")
135
  parser.add_argument("--output", type=Path, default=ROOT / "data/taxonomy.sqlite")
136
  parser.add_argument("--download", action="store_true", help="Download fresh NCBI metadata before rebuilding (about 2 GB).")
137
+ parser.add_argument("--coverage", type=Path, default=ROOT / "data/coverage.json", help="Published annotation inventory from refresh_coverage.py, used when present.")
138
  args = parser.parse_args()
139
  args.cache_dir.mkdir(parents=True, exist_ok=True)
140
  if args.download:
 
147
  out.write(chunk)
148
  temporary.replace(target)
149
  target.with_suffix(target.suffix + ".source.json").write_text(json.dumps(source))
150
+ build(args.cache_dir / "assembly_summary_genbank.txt", args.cache_dir / "taxdump.tar.gz", args.output, args.coverage)
151
 
152
 
153
  if __name__ == "__main__":
taxonomy.py CHANGED
@@ -15,6 +15,7 @@ class Taxonomy:
15
  self.path = Path(path)
16
  with closing(self.connect()) as conn:
17
  self.metadata = json.loads(conn.execute("SELECT value FROM metadata").fetchone()[0])
 
18
 
19
  def connect(self):
20
  conn = sqlite3.connect(f"{self.path.resolve().as_uri()}?mode=ro", uri=True)
@@ -36,16 +37,18 @@ class Taxonomy:
36
  root = conn.execute("SELECT * FROM taxa WHERE taxid=?", (int(taxid),)).fetchone()
37
  if root is None:
38
  raise ValueError("Taxon is absent from this assembly snapshot")
39
- ids, labels, parents, values, details = [], [], [], [], []
40
 
41
- def add(key, label, parent, count, rank):
42
  ids.append(str(key)); labels.append(label); parents.append(str(parent))
43
  values.append(count)
44
- details.append([rank, count / root['total_count'] * 100])
 
 
45
 
46
  def visit(row, parent="", depth=0):
47
  key = row['taxid']
48
- add(key, row['name'], parent, row['total_count'], row['rank'])
49
  if depth >= 3:
50
  return
51
  children = conn.execute("SELECT * FROM taxa WHERE parent_id=? AND taxid!=? ORDER BY total_count DESC", (key, key)).fetchall()
@@ -54,9 +57,10 @@ class Taxonomy:
54
  visit(child, key, depth + 1)
55
  other = sum(c['total_count'] for c in children[8:])
56
  if other:
57
- add(f"other-{key}", f"Other taxa ({len(children) - 8:,})", key, other, "Grouped for display; search to explore")
 
58
  if row['direct_count'] and children:
59
- add(f"direct-{key}", "Assigned directly to this taxon", key, row['direct_count'], "direct assignments")
60
 
61
  visit(root)
62
  children = conn.execute("SELECT * FROM taxa WHERE parent_id=? AND taxid!=? ORDER BY total_count DESC LIMIT 200", (root['taxid'], root['taxid'])).fetchall()
@@ -66,12 +70,21 @@ class Taxonomy:
66
  if current['taxid'] == 1:
67
  break
68
  current = conn.execute("SELECT * FROM taxa WHERE taxid=?", (current['parent_id'],)).fetchone()
 
 
 
 
 
 
 
69
  figure = go.Figure(go.Sunburst(ids=ids, labels=labels, parents=parents, values=values,
70
  branchvalues="total", customdata=details, sort=False,
71
- hovertemplate="<b>%{label}</b><br>%{customdata[0]}<br>%{value:,} assemblies<br>%{customdata[1]:.2f}% of selected group<extra></extra>"))
72
  figure.update_layout(height=620, margin=dict(t=10, b=10, l=10, r=10), template="plotly_white")
73
  count = root['total_count']
74
  note = f"### {root['name']}\n**{count:,} assemblies** · {count / self.metadata['assembly_count']:.2%} of this GenBank assembly snapshot\n\n" + " → ".join(reversed(lineage))
 
 
75
  choices = [(f"{r['name']} · {r['total_count']:,}", str(r['taxid'])) for r in children]
76
  return figure, note, choices, str(root['parent_id'])
77
 
@@ -84,7 +97,20 @@ def build_taxonomy_tab():
84
  meta = taxonomy.metadata
85
  gr.Markdown(f"## GenBank genome assemblies\n**{meta['assembly_count']:,} assemblies** · Snapshot built **{meta['created_at'][:10]}**\n\n"
86
  "Explore the NCBI taxonomy by assembly count. This inventory covers current GenBank genome assemblies; "
87
- "it does not count every nucleotide submission or indicate which genomes we have annotated.")
 
 
 
 
 
 
 
 
 
 
 
 
 
88
  initial, initial_note, initial_choices, _ = taxonomy.view("1")
89
  selected = gr.State("1")
90
  parent = gr.State("1")
@@ -102,7 +128,7 @@ def build_taxonomy_tab():
102
  gr.Markdown("Click chart sectors to expand the displayed branches; click the center to zoom out. "
103
  "Use **Explore a child group** or **Find a taxon** to load deeper branches. "
104
  "Each view includes three descendant levels and the eight largest children per node. "
105
- "**Other taxa** preserves the remaining assembly counts. Colors distinguish branches; they do not show annotation coverage.")
106
  with gr.Accordion("Snapshot sources and scope", open=False):
107
  gr.Markdown("Source: [NCBI GenBank assembly summary](https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt) "
108
  "and [NCBI Taxonomy](https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz). "
 
15
  self.path = Path(path)
16
  with closing(self.connect()) as conn:
17
  self.metadata = json.loads(conn.execute("SELECT value FROM metadata").fetchone()[0])
18
+ self.coverage = self.metadata.get("coverage")
19
 
20
  def connect(self):
21
  conn = sqlite3.connect(f"{self.path.resolve().as_uri()}?mode=ro", uri=True)
 
37
  root = conn.execute("SELECT * FROM taxa WHERE taxid=?", (int(taxid),)).fetchone()
38
  if root is None:
39
  raise ValueError("Taxon is absent from this assembly snapshot")
40
+ ids, labels, parents, values, details, colors = [], [], [], [], [], []
41
 
42
+ def add(key, label, parent, count, rank, covered):
43
  ids.append(str(key)); labels.append(label); parents.append(str(parent))
44
  values.append(count)
45
+ percent = 100 * covered / count if count else 0
46
+ colors.append(percent)
47
+ details.append([rank, count / root['total_count'] * 100, covered, percent])
48
 
49
  def visit(row, parent="", depth=0):
50
  key = row['taxid']
51
+ add(key, row['name'], parent, row['total_count'], row['rank'], row['covered_total'] if self.coverage else 0)
52
  if depth >= 3:
53
  return
54
  children = conn.execute("SELECT * FROM taxa WHERE parent_id=? AND taxid!=? ORDER BY total_count DESC", (key, key)).fetchall()
 
57
  visit(child, key, depth + 1)
58
  other = sum(c['total_count'] for c in children[8:])
59
  if other:
60
+ other_covered = sum(c['covered_total'] for c in children[8:]) if self.coverage else 0
61
+ add(f"other-{key}", f"Other taxa ({len(children) - 8:,})", key, other, "Grouped for display; search to explore", other_covered)
62
  if row['direct_count'] and children:
63
+ add(f"direct-{key}", "Assigned directly to this taxon", key, row['direct_count'], "direct assignments", row['covered_direct'] if self.coverage else 0)
64
 
65
  visit(root)
66
  children = conn.execute("SELECT * FROM taxa WHERE parent_id=? AND taxid!=? ORDER BY total_count DESC LIMIT 200", (root['taxid'], root['taxid'])).fetchall()
 
70
  if current['taxid'] == 1:
71
  break
72
  current = conn.execute("SELECT * FROM taxa WHERE taxid=?", (current['parent_id'],)).fetchone()
73
+ marker = dict(colors=colors, cmin=0, cmax=100, showscale=True,
74
+ colorscale=[[0, "#e2e8f0"], [0.000001, "#ccece6"], [0.1, "#80cdc1"], [0.5, "#279e91"], [1, "#005a50"]],
75
+ colorbar=dict(title="Assemblies with<br>annotations (%)", tickvals=[0, 25, 50, 75, 100]),
76
+ line=dict(color="white", width=1)) if self.coverage else {}
77
+ hover = "<b>%{label}</b><br>%{customdata[0]}<br>%{value:,} assemblies<br>%{customdata[1]:.2f}% of selected group"
78
+ if self.coverage:
79
+ hover += "<br><b>%{customdata[2]:,} with published annotations</b><br>%{customdata[3]:.2f}% of assemblies in this group"
80
  figure = go.Figure(go.Sunburst(ids=ids, labels=labels, parents=parents, values=values,
81
  branchvalues="total", customdata=details, sort=False,
82
+ marker=marker, hovertemplate=hover + "<extra></extra>"))
83
  figure.update_layout(height=620, margin=dict(t=10, b=10, l=10, r=10), template="plotly_white")
84
  count = root['total_count']
85
  note = f"### {root['name']}\n**{count:,} assemblies** · {count / self.metadata['assembly_count']:.2%} of this GenBank assembly snapshot\n\n" + " → ".join(reversed(lineage))
86
+ if self.coverage:
87
+ note += f"\n\n**{root['covered_total']:,} / {count:,} assemblies with published annotations ({root['covered_total'] / count:.2%})**. Assemblies may be partially annotated."
88
  choices = [(f"{r['name']} · {r['total_count']:,}", str(r['taxid'])) for r in children]
89
  return figure, note, choices, str(root['parent_id'])
90
 
 
97
  meta = taxonomy.metadata
98
  gr.Markdown(f"## GenBank genome assemblies\n**{meta['assembly_count']:,} assemblies** · Snapshot built **{meta['created_at'][:10]}**\n\n"
99
  "Explore the NCBI taxonomy by assembly count. This inventory covers current GenBank genome assemblies; "
100
+ "it does not count every nucleotide submission.")
101
+ if taxonomy.coverage:
102
+ coverage = taxonomy.coverage
103
+ gr.Markdown(f"**Annotation inventory updated {coverage['created_at'][:10]}** · "
104
+ f"**{coverage['matched_current_assemblies']:,} assemblies** have published annotations\n\n"
105
+ "**Area = assembly count. Color = percentage of assemblies with published annotations.** "
106
+ "Gray means no published annotations found in this inventory; darker teal means a higher fraction. "
107
+ "An assembly counts once if it has annotations, including partial assemblies. This is not the percentage of bases annotated. "
108
+ "The Annotation lookup tab currently searches a smaller subset.")
109
+ if coverage['unmatched_assembly_versions']:
110
+ gr.Markdown(f"{coverage['unmatched_assembly_versions']:,} published assembly versions are outside the current GenBank snapshot "
111
+ "and are excluded from the chart's percentages. Matching uses exact accession versions.")
112
+ else:
113
+ gr.Markdown("Annotation coverage inventory is unavailable; colors distinguish taxonomic branches only.")
114
  initial, initial_note, initial_choices, _ = taxonomy.view("1")
115
  selected = gr.State("1")
116
  parent = gr.State("1")
 
128
  gr.Markdown("Click chart sectors to expand the displayed branches; click the center to zoom out. "
129
  "Use **Explore a child group** or **Find a taxon** to load deeper branches. "
130
  "Each view includes three descendant levels and the eight largest children per node. "
131
+ "**Other taxa** preserves the remaining assembly counts and combines their annotation counts before calculating its percentage.")
132
  with gr.Accordion("Snapshot sources and scope", open=False):
133
  gr.Markdown("Source: [NCBI GenBank assembly summary](https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt) "
134
  "and [NCBI Taxonomy](https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz). "