Skip to content

Commit 15c193a

Browse files
jonn-smithjonn-smith
committed
Add analysis/ folder with hvp_viral_viz package
Introduce analysis/ as the home for downstream tools and analyses that operate on HVP-LR pipeline outputs but are not executed by the WDL workflow itself. First inhabitant is the hvp_viral_viz package — virus×ORF embedding, leiden clustering with rank-coherent ICTV resolution selection, biologically-labelled UMAPs, and per-cluster marker proteins from UniProt — consuming HvpViralProteinAnnotation outputs. - analysis/README.md documents conventions (one self-contained subdir per project, no bulk data in git, pinned deps, upstream-contract docs, smoke tests required). - Root README.md grows an analysis/ entry in the layout diagram and the "Where to put new things" table, plus a section linking to analysis/README.md. - hvp_viral_viz brings its own README, PLAN/HOST_FILL/THRESHOLD design docs, and all package modules (build_bfvd_refs, build_host_table, host_taxid_lookup, uniprot_host, uniprot_protein_name, ingest, plots, plots_tier2, scan_resolution, label_clusters, cluster_markers).
1 parent 33da9a6 commit 15c193a

18 files changed

Lines changed: 3549 additions & 0 deletions

‎README.md‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,9 @@ HVP-LR/
1313
│ ├── structs/Structs.wdl # Shared RuntimeAttr struct
1414
│ ├── tasks/<Category>/ # Task modules (Preprocessing, QC, Utility, ...)
1515
│ └── pipelines/<Tech>/<Cat>/ # End-to-end workflows (<Tech> = PacBio | ONT | ILMN | TechAgnostic)
16+
├── analysis/ # Downstream tools and analyses consuming
17+
│ └── <project>/ # pipeline outputs (figures, cluster labels,
18+
│ # marker tables, …). NOT executed by WDL.
1619
├── docs/ # Project docs (WDL style rules, etc.)
1720
├── scripts/ # CI helper scripts
1821
├── .github/workflows/ # docker.yml (build/push/scan), cd.yml (release)
@@ -36,6 +39,20 @@ HVP-LR/
3639
| A new pipeline / workflow | `wdl/pipelines/<Tech>/<Category>/<Name>.wdl` |
3740
| A workflow's inputs JSON | Alongside the WDL as `<Name>.inputs.json` |
3841
| Helper / one-off scripts | `scripts/` (avoid if not CI-related; otherwise put inside the relevant `docker/<image>/`) |
42+
| A downstream analysis / visualization tool | `analysis/<project>/` — see [analysis/README.md](analysis/README.md) for conventions |
43+
44+
## The `analysis/` folder
45+
46+
[`analysis/`](analysis/) holds downstream tools and analyses that operate on
47+
data produced by HVP-LR pipelines (and sibling reference datasets like
48+
BFVD / ICTV). It is **not** part of the WDL execution path — code here is
49+
invoked by humans (or notebooks) against pipeline outputs to produce
50+
figures, cluster labels, marker tables, and other interpretive artifacts.
51+
52+
Each project lives in its own self-contained subdirectory with its own
53+
`README.md`, dependency manifest, source, and tests. See
54+
[`analysis/README.md`](analysis/README.md) for full conventions and the
55+
list of current projects.
3956

4057
## Where to read next
4158

‎analysis/README.md‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
# analysis/
2+
3+
Tools and analyses that operate on data produced by the HVP-LR pipelines.
4+
5+
This folder is **not** part of the WDL workflow execution path. WDL tasks
6+
live under [`wdl/`](../wdl/) and container images under [`docker/`](../docker/).
7+
What lives here instead is downstream code that consumes pipeline outputs
8+
(hits TSVs, AnnData files, BFVD reference joins, foldseek m8 tables, etc.)
9+
and turns them into figures, reports, cluster labels, marker tables, and
10+
other interpretive artifacts.
11+
12+
## When to add something here
13+
14+
Put a new project under `analysis/<name>/` when **all** of these are true:
15+
16+
- It reads, transforms, plots, or annotates data emitted by an HVP-LR
17+
workflow (or a sibling reference dataset like BFVD / ICTV).
18+
- It does **not** need to run inside a Cromwell task as part of a sample's
19+
primary pipeline. If it needs to run on every sample as part of the
20+
pipeline, it belongs in `wdl/` with a container in `docker/`.
21+
- It is meaningful at the project level — useful to more than one sample,
22+
more than one user, or more than one experiment. One-off scratch
23+
notebooks do not belong here; check them into your own workspace.
24+
25+
## Conventions
26+
27+
Each project lives in its own subdirectory and is self-contained:
28+
29+
```
30+
analysis/<name>/
31+
├── README.md # what it does, how to run it, inputs/outputs
32+
├── pyproject.toml # if Python, for `pip install -e .`
33+
├── <package>/ # source (or src/<package>/ — either layout is fine)
34+
├── tests/ # pytest-style unit tests
35+
└── docs/ # any longer-form design docs
36+
```
37+
38+
Rules:
39+
40+
- **No bulk data in git.** Reference parquets, AnnData files, foldseek m8
41+
outputs, per-sample plots, and UniProt caches do **not** live here. Stash
42+
them in the relevant GCS bucket and document the URI in the project's
43+
README. Add `*.parquet`, `*.h5ad`, `out/`, `refs/cache/`, etc. to a local
44+
`.gitignore` inside the project.
45+
- **Pin your dependencies.** A `pyproject.toml` (or `requirements.txt` with
46+
hashes) keeps the analysis reproducible.
47+
- **Document the upstream contract.** State which workflow's outputs you
48+
consume and which columns / schema you assume. If the upstream schema
49+
changes, the analysis must update in the same PR.
50+
- **Tests are not optional.** At minimum a smoke test per module that
51+
loads a tiny fixture and runs the main entry point.
52+
53+
## Current projects
54+
55+
| Project | Path | Purpose |
56+
|---------|------|---------|
57+
| `hvp_viral_viz` | [`analysis/hvp_viral_viz/`](hvp_viral_viz/) | Downstream visualization of `HvpViralProteinAnnotation` outputs: virus×ORF embedding, leiden clustering with rank-coherent resolution selection, biologically-labelled UMAPs, per-cluster marker proteins from UniProt. Consumes `hits.tsv` + BFVD reference joins. |

‎analysis/hvp_viral_viz/HOST_FILL.md‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
# Host species fill — alternatives + chosen strategy
2+
3+
## Problem
4+
5+
`annotated_hits.tsv` carries no host info. BFVD reference tables
6+
(`bfvd/`) carry virus taxonomy but no host either. Need
7+
`taxid → host_group` mapping for UMAP-by-host plot.
8+
9+
## Definition of `host_group`
10+
11+
Coarse buckets that compress 23,614 viral taxids into ~10 categories
12+
useful for visual stratification:
13+
14+
| host_group | examples |
15+
|---|---|
16+
| `human` | HIV, HSV, HPV, SARS-CoV-2, common-cold coronaviruses |
17+
| `vertebrate_nonhuman` | bat coronaviruses, avian influenza, fish iridoviruses |
18+
| `arthropod` | insect baculoviruses, mosquito-vectored viruses (vector not host of replication) |
19+
| `plant` | TMV, plant geminiviruses |
20+
| `fungus` | mycoviruses |
21+
| `protist` | algal viruses, mimivirus host *Acanthamoeba* |
22+
| `bacteria` | bacteriophages (Caudoviricetes, Microviridae, etc.) |
23+
| `archaea` | archaeal viruses |
24+
| `metazoan_other` | nematodes, molluscs (abalone herpesvirus) |
25+
| `unknown` | uncharacterized environmental virus |
26+
27+
`host_organism` retained as free-text alongside (e.g. "Homo sapiens",
28+
"Gordonia bronchialis") for hover-text in interactive plots.
29+
30+
## Source ranking — best-first
31+
32+
Use a layered fill: each layer fills any taxid the previous left
33+
unknown, with `host_source` recording which layer landed each value.
34+
35+
### 1. ICTV Virus Metadata Resource (VMR) — primary
36+
37+
- Single authoritative table per ICTV release (MSL40 current as of 2026).
38+
- Columns include `Host source` (human-curated) at species / genus / family granularity.
39+
- Pull URL: `https://ictv.global/vmr/current` → CSV → ~17k rows mapping ICTV taxon → host source.
40+
- **Limit:** ICTV taxonomy ≠ NCBI taxid. Need ICTV name → NCBI taxid bridge.
41+
Two paths:
42+
- Use ICTV's `Genome composition` + `Virus name(s)` to fuzzy-match to NCBI taxonomy names (provided in BFVD lineage already).
43+
- Use `ictv-taxonomy-id` field in VMR (if released after MSL38, which carries NCBI taxid for ~80% of entries).
44+
- Expected coverage: 60-75% of BFVD taxids hit ICTV exactly at species or genus level.
45+
46+
### 2. UniProt REST `organism_host` — secondary
47+
48+
- Each viral UniProt entry has an `Organism host` cross-reference (when
49+
characterized). Query API at
50+
`https://rest.uniprot.org/uniprotkb/{accession}.json?fields=organism_host`.
51+
- Coverage on BFVD: estimated 30-50% (many BFVD entries are uncharacterized
52+
metagenome predictions where UniProt has no host annotation).
53+
- For the unique 347k uniprots in BFVD, this is ~70k REST calls if
54+
batched (UniProt allows 500 ids per stream/search query → ~700 calls
55+
with rate-limit pause). Cacheable as `refs/uniprot_host.parquet`.
56+
57+
### 3. NCBI Taxonomy host field — tertiary
58+
59+
- NCBI taxonomy entries carry `host` in the dataset (efetch from
60+
`taxonomy` db). Coverage similar to UniProt but for the taxid (not
61+
uniprot) — fills viruses where the species is known but the UniProt
62+
entry is bare.
63+
- ~24k taxids → 24k efetch calls (cacheable; batched in groups of 200).
64+
65+
### 4. Rule-based by lineage — quaternary
66+
67+
- Sweep remaining unknowns using BFVD lineage we already have.
68+
Mapping rules from clade → host_group:
69+
70+
```
71+
class Caudoviricetes → bacteria
72+
family Microviridae → bacteria
73+
family Inoviridae → bacteria
74+
family Tectiviridae → bacteria
75+
family Plasmaviridae → bacteria
76+
family Corticoviridae → bacteria
77+
family Sphaerolipoviridae → bacteria / archaea (mostly bacteria; flag dual)
78+
family Fuselloviridae → archaea
79+
family Rudiviridae → archaea
80+
family Lipothrixviridae → archaea
81+
family Globuloviridae → archaea
82+
family Ampullaviridae → archaea
83+
family Bicaudaviridae → archaea
84+
family Salterprovirus → archaea
85+
order Geplafuvirales → plant
86+
family Geminiviridae → plant
87+
family Caulimoviridae → plant
88+
family Bromoviridae → plant
89+
family Closteroviridae → plant
90+
family Tymoviridae → plant
91+
family Tobamoviridae → plant
92+
family Virgaviridae → plant
93+
family Mimiviridae → protist (Acanthamoeba)
94+
family Marseilleviridae → protist
95+
family Phycodnaviridae → protist (algae)
96+
family Pithoviridae → protist
97+
family Asfarviridae → vertebrate_nonhuman (mostly suid)
98+
family Iridoviridae → vertebrate_nonhuman + arthropod (flag mixed)
99+
family Baculoviridae → arthropod
100+
family Polydnaviridae → arthropod
101+
family Nudiviridae → arthropod
102+
family Hytrosaviridae → arthropod
103+
family Hytrosaviridae → arthropod
104+
family Partitiviridae → fungus / plant (flag mixed)
105+
family Totiviridae → fungus / protist (flag mixed)
106+
family Hypoviridae → fungus
107+
family Polyomaviridae → vertebrate_nonhuman (mostly mammals; human via taxid override)
108+
family Papillomaviridae → vertebrate_nonhuman (mostly mammals; human via taxid override)
109+
family Herpesviridae → vertebrate_nonhuman (human via taxid override)
110+
family Coronaviridae → vertebrate_nonhuman (human via taxid override)
111+
family Flaviviridae → vertebrate_nonhuman + arthropod (vector; flag)
112+
family Filoviridae → vertebrate_nonhuman (human via taxid override)
113+
family Retroviridae → vertebrate_nonhuman (human via taxid override)
114+
```
115+
116+
- After-rule sweep, force any taxid in NCBI `Viruses_human` curated
117+
list → `host_group = human` (overrides family default).
118+
119+
### 5. Scientific-name pattern match — last resort
120+
121+
- BFVD's scientific name often embeds the host: `Gordonia phage Dardanus`,
122+
`Abalone herpesvirus`, `Acanthamoeba polyphaga mimivirus`. Regex:
123+
- `phage` or `prophage` → bacteria
124+
- First word is a known prokaryote genus → bacteria (via NCBI taxonomy
125+
dump cross-check)
126+
- First word ends in `virus` and is preceded by a host species name →
127+
derive host_organism from preceding tokens
128+
- This catches an estimated 20-30% of remaining unknowns. Tagged
129+
`host_source = name_regex` so we can downweight in any quantitative
130+
analysis.
131+
132+
## Alternative sources NOT used (and why)
133+
134+
- **ViralZone (Expasy)** — hand-curated, web-only, no bulk download.
135+
Useful for spot-checks of high-interest viruses, not bulk fill.
136+
- **NCBI Virus** (the portal, not the taxdump) — overlaps with NCBI
137+
taxonomy host field; portal API is rate-limited and intended for
138+
interactive use.
139+
- **CHVD / IMG/VR / GOV2** — environmental virus catalogs with host
140+
prediction. Useful if our novel hits cluster around uncharacterized
141+
BFVD entries; would attach via target uniprot → environmental virus
142+
catalog ID → predicted host. Out of scope for v1.
143+
- **CRISPR-spacer-based host prediction** (e.g. iPHoP, PHIST) — strongest
144+
predictor for phages whose host has CRISPR records. Worth running on
145+
the BFVD subset that fails layers 1-5; deferred to a separate task
146+
since it requires downloading several large reference DBs.
147+
- **Sequence-similarity-based host transfer** (BLAST our query against
148+
a host-tagged virus DB, transfer host of nearest neighbor) — circular
149+
for our use case (we'd be using the same hits we're trying to label).
150+
151+
## Cost of layers 1-3
152+
153+
| Layer | API calls | Wall time (single-thread) | Cache reuse |
154+
|---|---|---|---|
155+
| 1 ICTV VMR | 1 (CSV download) | 5 s | Forever |
156+
| 2 UniProt REST | ~700 batched queries | ~30 min (with rate-limit) | Per uniprot |
157+
| 3 NCBI efetch | ~120 batches of 200 | ~10 min | Per taxid |
158+
159+
All caches stored under `refs/cache/` as parquet; only refetch on
160+
explicit `--refresh-cache`.
161+
162+
## Expected final coverage
163+
164+
| Source layer | Estimated cumulative coverage |
165+
|---|---|
166+
| ICTV VMR | 60-75% |
167+
| + UniProt host | 75-85% |
168+
| + NCBI taxonomy host | 80-88% |
169+
| + Rule-based by lineage | 92-97% |
170+
| + Name regex | 95-99% |
171+
| `unknown` residual | 1-5% |
172+
173+
All `unknown` entries plot as gray in the host UMAP and are excluded
174+
from host-stratified abundance bars (counted separately).
175+
176+
## Audit trail
177+
178+
`host_fill_audit.tsv` columns:
179+
```
180+
taxid, scientific_name, family, genus, host_group, host_organism,
181+
host_source, host_confidence
182+
```
183+
184+
`host_confidence` ∈ {`high`, `medium`, `low`}:
185+
- `high` = layers 1-3 (curated)
186+
- `medium` = layer 4 (rule by lineage)
187+
- `low` = layer 5 (name regex)
188+
189+
Plots in `out/` can be regenerated against `host_confidence ≥ medium`
190+
to drop name-regex entries when reviewer wants conservative figures.
191+
192+
## Chosen v1 scope
193+
194+
Implement layers 1, 4, 5 first (zero cost, fastest signal). Add
195+
layers 2 + 3 once layers 1+4+5 leave a residual we can quantify. This
196+
keeps `make refs` fast for iteration and avoids burning UniProt /
197+
NCBI quota during dev.

‎analysis/hvp_viral_viz/PLAN.md‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
# Implementation plan — viral-viz
2+
3+
Sequenced. Each phase produces a concrete artifact + check.
4+
5+
## Phase 0 — Foundations (done as part of this plan write-up)
6+
7+
- [x] Folder skeleton + docs (README, PLAN, THRESHOLD, HOST_FILL)
8+
- [x] Confirm BFVD has taxonomy lineage with realm→species (351,125 rows)
9+
- [x] Confirm host species NOT in BFVD; external fill needed
10+
- [x] Empirical sample of foldseek m8 score distribution
11+
- [x] Pick threshold + write rationale (see THRESHOLD.md)
12+
13+
## Phase 1 — Reference tables (refs/)
14+
15+
1. **`build_bfvd_refs.py`** — collapse BFVD into per-uniprot tables.
16+
Inputs: 3 BFVD TSVs from `/workspace/bfvd/`.
17+
Outputs in `refs/`:
18+
- `uniprot_taxid.parquet` — uniprot_id, taxid (≈ 347,514 rows)
19+
- `taxid_lineage.parquet` — taxid, scientific_name, realm, kingdom, phylum, class, order, family, genus, species (≈ 23,614 rows; parse `d_/k_/p_/c_/o_/f_/g_/s_` prefixes)
20+
- `uniprot_qc.parquet` — uniprot_id, model_id, avg_plddt, ptm, splitted, version (best per uniprot when split)
21+
Validate: every BFVD uniprot has taxid; ≥ 99% of taxids resolve to realm + species.
22+
23+
2. **`build_host_table.py`** — produce taxid→host_group.
24+
Strategy in `HOST_FILL.md`. Outputs:
25+
- `taxid_host.parquet` — taxid, host_group, host_source (rule|ictv|uniprot|ncbi|unknown), host_organism
26+
- `host_fill_audit.tsv` — per-taxid which source filled it
27+
28+
## Phase 2 — Hit ingestor (src/hvp_viral_viz/ingest.py)
29+
30+
3. **`ingest_run`** — CLI consuming one m8 + sample metadata:
31+
- Parse m8 (7 cols: query, target, evalue, bits, fident, alnlen, mismatch)
32+
- Apply discovery threshold (THRESHOLD.md → bits ≥ 50 AND evalue ≤ 1e-5)
33+
- Join `target` → uniprot_id → taxid → lineage → host
34+
- Tag each query with ORF source from query prefix (`vs2|`, `assembly|`, `rescued|`)
35+
- Emit `out/per_sample/<sample>.hits.parquet`
36+
37+
4. **`build_anndata`** — combine N per-sample parquets:
38+
- **Framing A (cohort, primary):** virus × sample matrix
39+
- obs = unique target uniprots (or collapsed to taxid for sparser viz)
40+
- var = sample_id (current = 1 sample)
41+
- .X = sum bits (or count hits, or binary) per (target, sample)
42+
- .obs columns: taxid, scientific_name, family, genus, host_group, avg_plddt, ptm
43+
- .var columns: sample_id, entity_type, cohort, donor, n_queries, n_orfs_vs2, n_orfs_assembly, n_orfs_rescued
44+
- **Framing B (single-sample fallback):** virus × query matrix
45+
- obs = unique targets, var = unique queries, .X = bits
46+
- Only computed for samples flagged as singleton viz target
47+
- Save as `out/anndata/cohort.h5ad` + `out/anndata/<sample>.singleton.h5ad`
48+
49+
## Phase 3 — Tier 1 plots (user-requested)
50+
51+
5. **Count matrix** — `plot_count_heatmap.py`
52+
- Clustered heatmap (viruses × samples), rows=top-N by total bits
53+
- CSV companion: `out/tables/virus_x_sample_counts.csv`
54+
55+
6. **UMAP × taxonomy** — `plot_umap_taxonomy.py`
56+
- Cohort mode (Framing A): `sc.pp.normalize` → `sc.pp.pca` → `sc.pp.neighbors` → `sc.tl.umap`; color by family / genus
57+
- Single mode (Framing B): same pipeline on virus×query matrix
58+
- One PNG per taxonomy rank: realm, phylum, class, family, genus
59+
60+
7. **UMAP × host** — `plot_umap_host.py`
61+
- Reuse coords from #6; recolor by host_group from `taxid_host`
62+
63+
8. **Abundance bar** — `plot_abundance_bar.py`
64+
- Top-N (default 30) viruses by hit count, stacked by family
65+
- Also: top-N by summed bits
66+
67+
## Phase 4 — Tier 2 diagnostic + interpretive plots
68+
69+
9. **Krona / sunburst** — `plot_taxonomy_sunburst.py` (plotly)
70+
10. **ORF-source Sankey** — `plot_orf_source_sankey.py` (vs2/assembly/rescued → cluster → family)
71+
11. **bits×fident scatter** — `plot_quality_scatter.py` colored by family
72+
12. **plDDT distribution per family** — `plot_plddt_per_family.py`
73+
13. **Hit-coverage histogram** — `plot_coverage_histogram.py` (alnlen / query_len)
74+
75+
## Phase 5 — Tier 3 biology-first (deferred unless asked)
76+
77+
14. Per-genome reconstruction grid — only when ≥ 1 target has many query hits
78+
15. Phylogenetic placement — out of scope for v1
79+
16. Bipartite network — out of scope for v1
80+
81+
## Phase 6 — CLI + reproducibility
82+
83+
17. **`hvp-viral-viz` console entry** — `pip install -e .` exposes:
84+
```
85+
hvp-viral-viz refs build
86+
hvp-viral-viz ingest <m8> --sample <id> --entity-type <type>
87+
hvp-viral-viz anndata build --out cohort.h5ad
88+
hvp-viral-viz plot all --anndata cohort.h5ad --out out/
89+
```
90+
18. **`Makefile`** at viral-viz root with `make refs`, `make ingest`, `make plots`, `make all`
91+
19. **`environment.yml`** — scanpy, anndata, plotly, pyarrow, requests
92+
20. **`tests/`** — unit tests on small synthetic m8 + small BFVD subset
93+
94+
## Phase 7 — Future cohort
95+
96+
21. Re-run `ingest_run` per new sample; `build_anndata` rebuilds cohort h5ad
97+
22. Sample-level PCoA + diversity plots gated on n_samples ≥ 3
98+
23. Per-virus prevalence vs mean-bits plot gated on n_samples ≥ 5
99+
100+
## Open questions to close before coding Phase 2
101+
102+
- best-hit per (query, target) or sum all (query, target) bits? → sum bits in `.X`, also store binary 1/0 layer
103+
- count "virus" at uniprot, taxid, or species level? → keep uniprot as primary obs key, derive species-level rollup as secondary AnnData
104+
- treat splitted models (`A0A2P1GMZ4_1` vs `A0A2P1GMZ4`) as same virus? → yes; collapse on uniprot stem

0 commit comments

Comments
 (0)