Add GenBank assembly taxonomy tab and refreshable SQLite snapshot
Browse files- .gitattributes +1 -0
- FEATURE_BACKLOG.md +5 -4
- README.md +19 -1
- app.py +42 -37
- data/taxonomy.sqlite +3 -0
- refresh_taxonomy.py +129 -0
- taxonomy.py +126 -0
.gitattributes
CHANGED
|
@@ -33,3 +33,4 @@ saved_model/**/* filter=lfs diff=lfs merge=lfs -text
|
|
| 33 |
*.zip filter=lfs diff=lfs merge=lfs -text
|
| 34 |
*.zst filter=lfs diff=lfs merge=lfs -text
|
| 35 |
*tfevents* filter=lfs diff=lfs merge=lfs -text
|
|
|
|
|
|
| 33 |
*.zip filter=lfs diff=lfs merge=lfs -text
|
| 34 |
*.zst filter=lfs diff=lfs merge=lfs -text
|
| 35 |
*tfevents* filter=lfs diff=lfs merge=lfs -text
|
| 36 |
+
data/taxonomy.sqlite filter=lfs diff=lfs merge=lfs -text
|
FEATURE_BACKLOG.md
CHANGED
|
@@ -9,7 +9,7 @@ Last updated: 2026-09-30.
|
|
| 9 |
| ID | Feature | Status |
|
| 10 |
| --- | --- | --- |
|
| 11 |
| F001 | GenBank inference coverage map | Planned; not implemented |
|
| 12 |
-
| F002 | Taxonomy map and navigation |
|
| 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 |
|
|
@@ -36,6 +36,8 @@ Decisions and dependencies for implementation:
|
|
| 36 |
|
| 37 |
## F002 — Taxonomy map and navigation
|
| 38 |
|
|
|
|
|
|
|
| 39 |
**Goal:** Provide a taxonomy-based view of the available genomes and annotations.
|
| 40 |
|
| 41 |
Desired behavior:
|
|
@@ -47,9 +49,8 @@ Desired behavior:
|
|
| 47 |
|
| 48 |
Decisions and dependencies for implementation:
|
| 49 |
|
| 50 |
-
-
|
| 51 |
-
-
|
| 52 |
-
- Define how to represent missing, unresolved, or changed taxon IDs.
|
| 53 |
|
| 54 |
## F003 — Comparison with RefSeq annotations
|
| 55 |
|
|
|
|
| 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 |
|
|
|
|
| 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 |
|
| 43 |
Desired behavior:
|
|
|
|
| 49 |
|
| 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 |
|
README.md
CHANGED
|
@@ -15,6 +15,24 @@ Search an indexed subset of [HuggingFaceBio/genbank-annotations](https://hugging
|
|
| 15 |
|
| 16 |
Future features and open design decisions are tracked in the [feature backlog](FEATURE_BACKLOG.md).
|
| 17 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 18 |
## Run locally
|
| 19 |
|
| 20 |
```bash
|
|
@@ -30,7 +48,7 @@ Open http://localhost:7860. Remote reads require a Hugging Face login (`hf auth
|
|
| 30 |
|
| 31 |
The preview Space is [HuggingFaceBio/genbank-annotation-explorer](https://huggingface.co/spaces/HuggingFaceBio/genbank-annotation-explorer), with private visibility.
|
| 32 |
|
| 33 |
-
To deploy or update it, upload
|
| 34 |
|
| 35 |
## Remote prototype scope
|
| 36 |
|
|
|
|
| 15 |
|
| 16 |
Future features and open design decisions are tracked in the [feature backlog](FEATURE_BACKLOG.md).
|
| 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 |
+
|
| 26 |
+
### Refresh the taxonomy snapshot
|
| 27 |
+
|
| 28 |
+
```bash
|
| 29 |
+
python refresh_taxonomy.py --download
|
| 30 |
+
```
|
| 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 |
|
| 38 |
```bash
|
|
|
|
| 48 |
|
| 49 |
The preview Space is [HuggingFaceBio/genbank-annotation-explorer](https://huggingface.co/spaces/HuggingFaceBio/genbank-annotation-explorer), with private visibility.
|
| 50 |
|
| 51 |
+
To deploy or update it, upload the app code and documentation, including `data/catalog.sqlite` and `data/taxonomy.sqlite`; exclude `.cache/` and local environments. Set `HF_TOKEN` as a Space secret with read access to the bucket. Include `data/sample.parquet` and `data/manifest.json` if offline sample mode is desired. Credentials stay on the server and must never be committed. The README metadata configures the runtime, following the [Gradio Spaces documentation](https://huggingface.co/docs/hub/spaces-sdks-gradio). Keep the Space private for the initial preview.
|
| 52 |
|
| 53 |
## Remote prototype scope
|
| 54 |
|
app.py
CHANGED
|
@@ -13,6 +13,7 @@ import gradio as gr
|
|
| 13 |
import pyarrow.parquet as pq
|
| 14 |
import plotly.graph_objects as go
|
| 15 |
|
|
|
|
| 16 |
from catalog import Catalog
|
| 17 |
from remote_catalog import RemoteCatalog, RemoteReadError
|
| 18 |
|
|
@@ -119,43 +120,47 @@ def build_app(catalog=None):
|
|
| 119 |
|
| 120 |
with gr.Blocks(title="GenBank Annotation Explorer", delete_cache=(3600, 3600)) as demo:
|
| 121 |
gr.Markdown("# GenBank Annotation Explorer\nExplore predicted coding regions by **assembly or contig accession**.")
|
| 122 |
-
gr.
|
| 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 |
-
|
| 156 |
-
|
| 157 |
-
|
| 158 |
-
|
|
|
|
|
|
|
|
|
|
|
|
|
| 159 |
outputs = [status, results, segment, file]
|
| 160 |
for event in (search_button.click, accession.submit):
|
| 161 |
event(search, accession, outputs).then(select_segment, [segment, mode, threshold], [metadata, start, end, plot, note, file])
|
|
|
|
| 13 |
import pyarrow.parquet as pq
|
| 14 |
import plotly.graph_objects as go
|
| 15 |
|
| 16 |
+
from taxonomy import build_taxonomy_tab
|
| 17 |
from catalog import Catalog
|
| 18 |
from remote_catalog import RemoteCatalog, RemoteReadError
|
| 19 |
|
|
|
|
| 120 |
|
| 121 |
with gr.Blocks(title="GenBank Annotation Explorer", delete_cache=(3600, 3600)) as demo:
|
| 122 |
gr.Markdown("# GenBank Annotation Explorer\nExplore predicted coding regions by **assembly or contig accession**.")
|
| 123 |
+
with gr.Tabs():
|
| 124 |
+
with gr.Tab("GenBank taxonomy"):
|
| 125 |
+
build_taxonomy_tab()
|
| 126 |
+
with gr.Tab("Annotation lookup"):
|
| 127 |
+
gr.Markdown(f"**Prototype {scope} · {len(catalog.records):,} segments · "
|
| 128 |
+
f"{len(catalog.manifest['assemblies']):,} assemblies · {catalog.manifest['bases']:,} bases**\n\n"
|
| 129 |
+
"Source: [HuggingFaceBio/genbank-annotations](https://huggingface.co/buckets/HuggingFaceBio/genbank-annotations). "
|
| 130 |
+
f"These are model-predicted CDS probabilities, not curated gene features. Search covers the {scope} only. "
|
| 131 |
+
+ (f"Annotations load on demand from **{len(catalog.manifest['sources'])} bucket files**. "
|
| 132 |
+
f"Index updated {catalog.manifest['created_at'][:10]}." if remote_mode else ""))
|
| 133 |
+
with gr.Row():
|
| 134 |
+
accession = gr.Textbox(label="Accession ID", placeholder="Assembly (GCA_…) or contig accession", scale=5)
|
| 135 |
+
search_button = gr.Button("Find annotations", variant="primary", scale=1)
|
| 136 |
+
examples = [first["assembly_accession"], first["record_name"]]
|
| 137 |
+
if remote_mode:
|
| 138 |
+
with closing(catalog.connect()) as conn:
|
| 139 |
+
examples += [r[0] for r in conn.execute("SELECT record_name FROM segments WHERE id IN (SELECT min(id) FROM segments GROUP BY object_path)")]
|
| 140 |
+
examples += ["JBPJTW010000350.1"]
|
| 141 |
+
gr.Examples(examples=[[e] for e in dict.fromkeys(examples)], inputs=accession)
|
| 142 |
+
status = gr.Markdown("Enter an accession or select an example. IDs are case-insensitive; version suffixes are optional.")
|
| 143 |
+
results = gr.Dataframe(value=catalog.table([]), interactive=False, label="Matching segments")
|
| 144 |
+
segment = gr.Dropdown(choices=[], label="Segment to explore", interactive=True)
|
| 145 |
+
with gr.Row():
|
| 146 |
+
start = gr.Number(label="Start (0-based, inclusive)", precision=0)
|
| 147 |
+
end = gr.Number(label="End (exclusive)", precision=0)
|
| 148 |
+
view = gr.Button("Update region")
|
| 149 |
+
with gr.Row():
|
| 150 |
+
mode = gr.Radio(["Probabilities", "Binary labels"], value="Probabilities", label="Viewer mode")
|
| 151 |
+
threshold = gr.Slider(0, 1, value=0.5, step=0.01, label="CDS threshold (max strand P > threshold)")
|
| 152 |
+
plot = gr.Plot(label="CDS tracks")
|
| 153 |
+
note = gr.Markdown("Choose a segment to see its probability tracks.")
|
| 154 |
+
with gr.Accordion("Segment metadata and provenance", open=False):
|
| 155 |
+
metadata = gr.JSON(label="Source metadata")
|
| 156 |
+
with gr.Row():
|
| 157 |
+
download = gr.Button("Prepare segment download")
|
| 158 |
+
file = gr.File(label="Original segment annotations (Parquet)", interactive=False)
|
| 159 |
+
with gr.Accordion("Browse indexed accessions and coverage", open=False):
|
| 160 |
+
gr.Markdown(f"Showing the first {len(all_ids):,} indexed segments. Search an accession to find other indexed records. "
|
| 161 |
+
"An indexed file does not imply complete coverage of its assembly.")
|
| 162 |
+
gr.Dataframe(value=catalog.table(all_ids), interactive=False)
|
| 163 |
+
gr.JSON(value=catalog.manifest, label="Index provenance" if remote_mode else "Sample provenance")
|
| 164 |
outputs = [status, results, segment, file]
|
| 165 |
for event in (search_button.click, accession.submit):
|
| 166 |
event(search, accession, outputs).then(select_segment, [segment, mode, threshold], [metadata, start, end, plot, note, file])
|
data/taxonomy.sqlite
ADDED
|
@@ -0,0 +1,3 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
version https://git-lfs.github.com/spec/v1
|
| 2 |
+
oid sha256:69388d96ee05e254b8eb8530a62a71d3d9e4f89e397427c1b7f99c140d488b55
|
| 3 |
+
size 24084480
|
refresh_taxonomy.py
ADDED
|
@@ -0,0 +1,129 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
"""Build a small, atomic SQLite taxonomy snapshot from NCBI metadata."""
|
| 2 |
+
import argparse
|
| 3 |
+
from collections import Counter
|
| 4 |
+
from datetime import datetime, timezone
|
| 5 |
+
import hashlib
|
| 6 |
+
import json
|
| 7 |
+
from pathlib import Path
|
| 8 |
+
import sqlite3
|
| 9 |
+
import tarfile
|
| 10 |
+
import urllib.request
|
| 11 |
+
|
| 12 |
+
ROOT = Path(__file__).resolve().parent
|
| 13 |
+
SOURCES = {
|
| 14 |
+
"assembly_summary_genbank.txt": "https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt",
|
| 15 |
+
"taxdump.tar.gz": "https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz",
|
| 16 |
+
}
|
| 17 |
+
|
| 18 |
+
|
| 19 |
+
def dump_rows(archive, name):
|
| 20 |
+
with archive.extractfile(name) as stream:
|
| 21 |
+
for line in stream:
|
| 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:
|
| 30 |
+
if line.startswith("##"):
|
| 31 |
+
continue
|
| 32 |
+
if line.startswith("#"):
|
| 33 |
+
header = line.lstrip("# ").rstrip("\n").split("\t")
|
| 34 |
+
columns = {name: header.index(name) for name in ("assembly_accession", "taxid", "version_status")}
|
| 35 |
+
continue
|
| 36 |
+
row = line.rstrip("\n").split("\t")
|
| 37 |
+
if not row[columns["assembly_accession"]].startswith("GCA_") or row[columns["version_status"]] != "latest":
|
| 38 |
+
skipped += 1
|
| 39 |
+
continue
|
| 40 |
+
try:
|
| 41 |
+
taxid = int(row[columns["taxid"]])
|
| 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:
|
| 48 |
+
for row in dump_rows(archive, "nodes.dmp"):
|
| 49 |
+
taxid = int(row[0])
|
| 50 |
+
parents[taxid], ranks[taxid] = int(row[1]), row[2]
|
| 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]
|
| 72 |
+
names = {-1: "Unresolved taxonomy", 1: "GenBank genome assemblies"}
|
| 73 |
+
for row in dump_rows(archive, "names.dmp"):
|
| 74 |
+
taxid = int(row[0])
|
| 75 |
+
if taxid in totals and taxid != 1 and row[3] == "scientific name":
|
| 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()
|
| 82 |
+
with open(path, "rb") as stream:
|
| 83 |
+
for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""):
|
| 84 |
+
digest.update(chunk)
|
| 85 |
+
source = {"url": SOURCES[path.name], "sha256": digest.hexdigest(), "bytes": path.stat().st_size}
|
| 86 |
+
sidecar = path.with_suffix(path.suffix + ".source.json")
|
| 87 |
+
if sidecar.exists():
|
| 88 |
+
source.update(json.loads(sidecar.read_text()))
|
| 89 |
+
provenance[path.name] = source
|
| 90 |
+
metadata = {"created_at": datetime.now(timezone.utc).isoformat(), "sources": provenance,
|
| 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():
|
| 109 |
+
parser = argparse.ArgumentParser(description=__doc__)
|
| 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:
|
| 116 |
+
for filename, url in SOURCES.items():
|
| 117 |
+
target = args.cache_dir / filename
|
| 118 |
+
temporary = target.with_suffix(target.suffix + ".part")
|
| 119 |
+
with urllib.request.urlopen(url, timeout=120) as response, open(temporary, "wb") as out:
|
| 120 |
+
source = {"downloaded_at": datetime.now(timezone.utc).isoformat(), "last_modified": response.headers.get("Last-Modified")}
|
| 121 |
+
while chunk := response.read(8 * 1024 * 1024):
|
| 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__":
|
| 129 |
+
main()
|
taxonomy.py
ADDED
|
@@ -0,0 +1,126 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
"""Read-only taxonomy navigation, independent of the annotation lookup index."""
|
| 2 |
+
from contextlib import closing
|
| 3 |
+
import json
|
| 4 |
+
from pathlib import Path
|
| 5 |
+
import sqlite3
|
| 6 |
+
|
| 7 |
+
import gradio as gr
|
| 8 |
+
import plotly.graph_objects as go
|
| 9 |
+
|
| 10 |
+
DATABASE = Path(__file__).resolve().parent / "data/taxonomy.sqlite"
|
| 11 |
+
|
| 12 |
+
|
| 13 |
+
class Taxonomy:
|
| 14 |
+
def __init__(self, path=DATABASE):
|
| 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)
|
| 21 |
+
conn.row_factory = sqlite3.Row
|
| 22 |
+
return conn
|
| 23 |
+
|
| 24 |
+
def search(self, query):
|
| 25 |
+
query = str(query or "").strip()
|
| 26 |
+
with closing(self.connect()) as conn:
|
| 27 |
+
if query.lstrip("-").isdigit():
|
| 28 |
+
rows = conn.execute("SELECT * FROM taxa WHERE taxid=?", (int(query),)).fetchall()
|
| 29 |
+
else:
|
| 30 |
+
pattern = "%" + query.replace("\\", "\\\\").replace("%", "\\%").replace("_", "\\_") + "%"
|
| 31 |
+
rows = conn.execute("SELECT * FROM taxa WHERE name LIKE ? ESCAPE '\\' ORDER BY total_count DESC LIMIT 100", (pattern,)).fetchall()
|
| 32 |
+
return [(f"{r['name']} · {r['total_count']:,} assemblies · taxid {r['taxid']}", str(r['taxid'])) for r in rows]
|
| 33 |
+
|
| 34 |
+
def view(self, taxid):
|
| 35 |
+
with closing(self.connect()) as conn:
|
| 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()
|
| 52 |
+
# At most 8 children per node; exact residual counts stay visible.
|
| 53 |
+
for child in children[:8]:
|
| 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()
|
| 63 |
+
lineage, current = [], root
|
| 64 |
+
while True:
|
| 65 |
+
lineage.append(current['name'])
|
| 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 |
+
|
| 78 |
+
|
| 79 |
+
def build_taxonomy_tab():
|
| 80 |
+
if not DATABASE.exists():
|
| 81 |
+
gr.Markdown("Taxonomy snapshot unavailable. Run `python refresh_taxonomy.py --download` to create it.")
|
| 82 |
+
return
|
| 83 |
+
taxonomy = Taxonomy()
|
| 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")
|
| 91 |
+
with gr.Row():
|
| 92 |
+
home = gr.Button("All GenBank assemblies")
|
| 93 |
+
up = gr.Button("Parent group")
|
| 94 |
+
child = gr.Dropdown(choices=initial_choices, label="Explore a child group (largest 200)", interactive=True)
|
| 95 |
+
with gr.Row():
|
| 96 |
+
query = gr.Textbox(label="Find a taxon", placeholder="Scientific name or NCBI taxon ID")
|
| 97 |
+
find = gr.Button("Search taxonomy")
|
| 98 |
+
matches = gr.Dropdown(choices=[], label="Taxon search results (up to 100, largest first)", interactive=True)
|
| 99 |
+
search_status = gr.Markdown()
|
| 100 |
+
note = gr.Markdown(initial_note)
|
| 101 |
+
chart = gr.Plot(value=initial, label="GenBank taxonomy")
|
| 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). "
|
| 109 |
+
"Counts roll up each assembly's taxon through its ancestors. Merged taxon IDs are resolved; unresolved IDs remain in an explicit group.")
|
| 110 |
+
gr.JSON(meta)
|
| 111 |
+
|
| 112 |
+
def taxonomy_view(taxid):
|
| 113 |
+
figure, text, choices, parent_id = taxonomy.view(taxid or "1")
|
| 114 |
+
return figure, text, gr.Dropdown(choices=choices, value=None), str(taxid or "1"), parent_id
|
| 115 |
+
|
| 116 |
+
def taxonomy_search(text):
|
| 117 |
+
choices = taxonomy.search(text) if str(text or "").strip() else []
|
| 118 |
+
return gr.Dropdown(choices=choices, value=None), f"{len(choices)} matching taxa shown." if choices else "No matching taxa. Try a scientific name or taxon ID."
|
| 119 |
+
|
| 120 |
+
outputs = [chart, note, child, selected, parent]
|
| 121 |
+
child.input(taxonomy_view, child, outputs)
|
| 122 |
+
matches.input(taxonomy_view, matches, outputs)
|
| 123 |
+
home.click(lambda: taxonomy_view("1"), outputs=outputs)
|
| 124 |
+
up.click(taxonomy_view, parent, outputs)
|
| 125 |
+
find.click(taxonomy_search, query, [matches, search_status])
|
| 126 |
+
query.submit(taxonomy_search, query, [matches, search_status])
|