Skip to content

Commit cf76698

Browse files
edkerkedkerk
committed
Add KEGG artefact-build scripts and HMM-cutoff calibration docs
1 parent d9f100c commit cf76698

5 files changed

Lines changed: 673 additions & 0 deletions

File tree

‎docs/kegg_data_format.md‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
# KEGG relational-table storage format
2+
3+
This note records *why* raven-python stores its KEGG-derived relational tables as
4+
**gzipped TSV**, and what other options we deliberately deferred. It applies to
5+
the maintainer-built KEGG artefacts described in PLAN.md §2.3b — the `ko_reaction`,
6+
`organism_gene_ko`, KO-name, and reaction-flag tables.
7+
8+
The reference GEM itself is stored as **gzipped RAVEN/cobra YAML**
9+
(`reference_model.yml.gz`) — RAVEN-native and MATLAB-readable, gzipped to match the
10+
tables (the YAML I/O transparently gzips on a `.gz` suffix). On the real KEGG dump
11+
this is ~1.1 MB (vs ~30 MB as SBML) for the full 12k-reaction gene-free model.
12+
13+
End users do not build any of this: the published artefacts are fetched and cached
14+
under `~/.cache/raven-python/data/kegg-<version>/` by `ensure_data` (see
15+
`raven_python.data`), mirroring how binaries are provisioned.
16+
17+
## Decision (current)
18+
19+
- **Small tables** (`ko_reaction`, `ko_names`, `rxn_flags`): **gzipped TSV
20+
(`.tsv.gz`)**. Each is well under 1 MB, so compression choice is irrelevant;
21+
gzip keeps them MATLAB-native and dependency-free.
22+
- **The large `organism_gene_ko` table**: **xz-compressed TSV
23+
(`organism_gene_ko.tsv.xz`), with rows sorted by `(organism, gene)`**.
24+
25+
Why the large table differs. It carries KEGG's ~9M gene↔KO associations and
26+
dominates the artefact set (≈78 MB as unsorted gzipped TSV). Two cheap,
27+
stdlib-only changes cut that to ≈27 MB (2.9×):
28+
29+
1. **Sort by `(organism, gene)`** before writing. Gene IDs from one organism
30+
share long common prefixes (locus tags, numeric runs); sorting makes them
31+
adjacent so the compressor can fold them. This alone takes 78 → 48 MB and
32+
happens to match the by-organism query pattern in
33+
`get_kegg_model_for_organism`. The sort is an external merge sort bounded to
34+
`chunk_rows` in memory (see `stream_organism_gene_ko`), so it stays scalable.
35+
2. **xz instead of gzip** (Python stdlib `lzma`). Its larger dictionary captures
36+
cross-row redundancy gzip's 32 KB window misses: sorted + xz reaches ≈27 MB.
37+
38+
- **pandas reads/writes both with zero extra dependencies** — compression is
39+
inferred from the `.gz`/`.xz` suffix; `lzma` and `gzip` are both stdlib, so
40+
this works natively on Windows, macOS, and Linux with no external binary.
41+
- **MATLAB caveat:** `readtable` reads gzipped TSV after a `gunzip`, but MATLAB
42+
has no built-in xz decompressor. The small tables stay MATLAB-native; the
43+
large table needs an external `unxz` (or Java/`7-Zip`) before `readtable` on
44+
the MATLAB side. The xz file is raven-python's (Python) primary artefact; this
45+
trades a little MATLAB convenience on the one big file for a ~3× size cut.
46+
47+
## Options considered
48+
49+
| Format | Python cost | MATLAB cost | Notes |
50+
| --- | --- | --- | --- |
51+
| **Gzipped TSV** ✅ | none (stdlib/pandas) | none (`readtable`) | Universal, text, types re-specified on read. Chosen. |
52+
| Parquet | `pyarrow` or `fastparquet` (~40–60 MB wheel) as a `raven-python[kegg]` extra | needs ≥ R2019a (`parquetread`, native) | Smaller, faster, typed, columnar. Win mainly at scale / repeated random access. |
53+
| SQLite | none (stdlib `sqlite3`) | **needs Database Toolbox** | Rejected: the MATLAB-side toolbox requirement breaks the "same files, both languages, no extra deps" goal. |
54+
55+
## When to revisit
56+
57+
Reconsider Parquet (or SQLite) if any of these become true:
58+
59+
- The `organism_gene_ko` table grows large enough that load *time* (not just
60+
size — the sort+xz change above already addresses on-disk size) becomes a real
61+
bottleneck. The remaining inefficiency is that building one species' model
62+
still loads all ~9M rows; sorted order makes a `searchsorted`/row-group
63+
by-organism read the natural next step before reaching for Parquet.
64+
- We start doing repeated random-access / columnar reads rather than a single
65+
load-once-per-run pattern.
66+
- A typed, self-describing schema becomes valuable (TSV loses dtypes; they are
67+
re-specified on read).
68+
69+
If revisited, prefer **Parquet** over SQLite (no MATLAB toolbox dependency; MATLAB
70+
reads Parquet natively from R2019a). It could be offered as an optional
71+
`raven-python[kegg]` extra (pyarrow) alongside the TSV default, rather than replacing
72+
it — keeping the dependency-free path intact for users who don't opt in.

‎docs/kegg_hmm_cutoff_calibration.md‎

Lines changed: 203 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,203 @@
1+
# KEGG HMM-query cut-off calibration
2+
3+
This note records the measurements behind the default KO-assignment parameters in
4+
`reconstruction/kegg/query.py` (`assign_kos` / `get_kegg_model_from_sequences`,
5+
pipeline step 3b.5) and IMPROVEMENTS **K15**. It is the evidence for moving away
6+
from RAVEN's `1e-50` cut-off.
7+
8+
## What the parameters do
9+
10+
`assign_kos` turns an `hmmscan` KO×gene E-value matrix into gene→KO assignments
11+
in three steps:
12+
13+
1. **`cutoff`** — keep hits with `evalue <= cutoff`.
14+
2. **`min_score_ratio_ko`** — within a KO, drop genes whose
15+
`log(evalue)/log(best_evalue_in_KO) < min_score_ratio_ko`.
16+
3. **`min_score_ratio_g`** — within a gene, drop KOs whose
17+
`log(evalue)/log(best_evalue_for_gene) < min_score_ratio_g`.
18+
19+
## Method
20+
21+
- **Data:** KEGG release 118. Libraries: the maintainer-built `prokaryotes.hmm`
22+
(831 MB) and `eukaryotes.hmm` (692 MB), 90 %-clustered, FFT-NS-2/PartTree (K12).
23+
- **Queries:** each organism's full proteome, extracted from `genes.pep`.
24+
- **Ground truth:** the organism's *real* KEGG gene→KO links, from the
25+
`organism_gene_ko` table (restricted, as the table is, to reaction-linked KOs).
26+
- **Prediction:** `assign_kos` output, with the `organism:` prefix stripped from
27+
query gene IDs so they match the bare gene IDs in the ground truth.
28+
- **Metrics (gene→KO level):** precision = |pred ∩ truth| / |pred|,
29+
recall = |pred ∩ truth| / |truth|, F1. Reaction-level: `rxn_rec` = fraction of
30+
the organism's true reactions recovered (KO→reaction via `ko_reaction`);
31+
`rxn_novel` = predicted reactions **not** in the annotation set.
32+
- Reproduce with [`scripts/analyze_hmm_cutoffs.py`](../scripts/analyze_hmm_cutoffs.py).
33+
34+
### Important caveat
35+
36+
All four organisms are **present in the libraries' training set**, so their own
37+
sequences hit their KO profiles strongly and recall is an upper bound. The
38+
calibration is therefore *relative* (how the parameters trade off, and where
39+
RAVEN's default sits relative to the signal), not an absolute accuracy estimate.
40+
A genome genuinely absent from KEGG would be the next validation. Also note that
41+
`rxn_novel` / "precision < 1" partly reflects **legitimate homology** KEGG never
42+
annotated for that organism (paralogs, un-curated genes), not pure error — so the
43+
precision figures are a lower bound on real precision.
44+
45+
## Organisms
46+
47+
| code | organism | library | proteome (seqs) | true gene→KO pairs | true reactions |
48+
|---|---|---|---|---|---|
49+
| `sce` | *Saccharomyces cerevisiae* (budding yeast) | euk | 6021 | 841 | 1217 |
50+
| `cme` | *Cyanidioschyzon merolae* (red alga) | euk | 5010 | 709 | 1157 |
51+
| `eco` | *Escherichia coli* K-12 MG1655 | prok | 4288 | 1071 | 1548 |
52+
| `mge` | *Mycoplasmoides genitalium* G37 (minimal genome) | prok | 476 | 110 | 211 |
53+
54+
`sce`/`eco` are model organisms; `cme`/`mge` are lesser-studied, `mge`
55+
additionally being a small, divergent genome.
56+
57+
## 1. E-value separation (the key result)
58+
59+
`log10(E-value)` percentiles of the best hit per (gene, KO) pair, split by whether
60+
the pair is in the organism's annotation (**matched**) or not (**novel**). Smaller
61+
(more negative) = stronger hit.
62+
63+
| organism | group | n | p50 | p90 | p95 | p99 |
64+
|---|---|---|---|---|---|---|
65+
| `sce` | matched | 835 | −155 | −75 | −59 | −33 |
66+
| `sce` | novel | 9467 | −8 | −2 | −0 | 1 |
67+
| `cme` | matched | 704 | −133 | −63 | −47 | −25 |
68+
| `cme` | novel | 10170 | −8 | −2 | −2 | 0 |
69+
| `eco` | matched | 1070 | −142 | −69 | −57 | −36 |
70+
| `eco` | novel | 27357 | −7 | −2 | −1 | −0 |
71+
| `mge` | matched | 110 | −100 | −42 | −35 | −15 |
72+
| `mge` | novel | 1904 | −4 | −2 | −1 | −0 |
73+
74+
**Reading:** matched pairs cluster at E ≈ 1e-100…1e-155; even their weakest 1 %
75+
sit at 1e-15…1e-36. Novel pairs cluster at ≈1e-8. The two are separated by ~20
76+
orders of magnitude. RAVEN's **`1e-50` lies inside the *matched* tail** (between
77+
the matched p90 and p95 for most organisms; past p90 for `mge`), so it discards
78+
real-but-weakly-scoring annotations while gaining nothing against the (far weaker)
79+
noise.
80+
81+
## 2. Cut-off sweep
82+
83+
`min_score_ratio_ko = 0.3`, `min_score_ratio_g = 0.8` fixed; gene→KO precision /
84+
recall / F1 and reaction recovery vs the annotation.
85+
86+
### `sce`
87+
| cutoff | gKO prec | gKO rec | gKO F1 | rxn rec | rxn novel |
88+
|---|---|---|---|---|---|
89+
| 1e-10 | 0.57 | 0.98 | 0.72 | 0.99 | 334 |
90+
| 1e-20 | 0.65 | 0.98 | 0.78 | 0.97 | 283 |
91+
| 1e-30 | 0.72 | 0.97 | 0.83 | 0.97 | 216 |
92+
| 1e-50 | 0.78 | 0.95 | 0.86 | 0.96 | 157 |
93+
| 1e-70 | 0.81 | 0.91 | 0.86 | 0.91 | 113 |
94+
| 1e-100 | 0.84 | 0.76 | 0.80 | 0.79 | 68 |
95+
96+
### `cme`
97+
| cutoff | gKO prec | gKO rec | gKO F1 | rxn rec | rxn novel |
98+
|---|---|---|---|---|---|
99+
| 1e-10 | 0.50 | 0.98 | 0.67 | 1.00 | 541 |
100+
| 1e-20 | 0.57 | 0.98 | 0.72 | 1.00 | 421 |
101+
| 1e-30 | 0.61 | 0.97 | 0.75 | 0.97 | 367 |
102+
| 1e-50 | 0.70 | 0.93 | 0.80 | 0.94 | 307 |
103+
| 1e-70 | 0.75 | 0.85 | 0.80 | 0.87 | 223 |
104+
| 1e-100 | 0.80 | 0.70 | 0.75 | 0.71 | 136 |
105+
106+
### `eco`
107+
| cutoff | gKO prec | gKO rec | gKO F1 | rxn rec | rxn novel |
108+
|---|---|---|---|---|---|
109+
| 1e-10 | 0.53 | 0.99 | 0.69 | 0.99 | 382 |
110+
| 1e-20 | 0.57 | 0.99 | 0.73 | 0.99 | 300 |
111+
| 1e-30 | 0.60 | 0.98 | 0.75 | 0.99 | 268 |
112+
| 1e-50 | 0.67 | 0.95 | 0.78 | 0.98 | 198 |
113+
| 1e-70 | 0.73 | 0.88 | 0.80 | 0.93 | 157 |
114+
| 1e-100 | 0.82 | 0.74 | 0.77 | 0.80 | 96 |
115+
116+
### `mge`
117+
| cutoff | gKO prec | gKO rec | gKO F1 | rxn rec | rxn novel |
118+
|---|---|---|---|---|---|
119+
| 1e-10 | 0.52 | 0.98 | 0.68 | 0.99 | 75 |
120+
| 1e-20 | 0.62 | 0.96 | 0.75 | 0.98 | 51 |
121+
| 1e-30 | 0.65 | 0.95 | 0.77 | 0.98 | 39 |
122+
| 1e-50 | 0.77 | 0.84 | 0.80 | 0.87 | 29 |
123+
| 1e-70 | 0.85 | 0.73 | 0.78 | 0.73 | 27 |
124+
| 1e-100 | 0.87 | 0.50 | 0.64 | 0.47 | 21 |
125+
126+
**Reading:** recall is flat-and-high from 1e-10 to ~1e-30, then falls as the
127+
cut-off eats into the matched tail — gently for model organisms, sharply for the
128+
divergent `mge` (rxn recall 0.98 → 0.87 from 1e-30 → 1e-50, → 0.47 at 1e-100).
129+
The recall lost to a stricter cut-off is *not* noise rejection (noise is at 1e-8);
130+
it is real annotation. `rxn_novel` shrinks with stricter cut-offs because strong
131+
un-annotated homologs are also removed.
132+
133+
## 3. Score-ratio sweep (`cutoff = 1e-50`)
134+
135+
| organism | ko ratio | g ratio | gKO prec | gKO rec | gKO F1 |
136+
|---|---|---|---|---|---|
137+
| `sce` | 0.0 | 0.50 | 0.61 | 0.96 | 0.74 |
138+
| `sce` | 0.0 | 0.80 | 0.77 | 0.95 | 0.85 |
139+
| `sce` | 0.0 | 0.95 | 0.84 | 0.93 | 0.88 |
140+
| `sce` | 0.3 | 0.80 | 0.78 | 0.95 | 0.86 |
141+
| `sce` | 0.5 | 0.80 | 0.80 | 0.95 | 0.86 |
142+
| `cme` | 0.0 | 0.50 | 0.53 | 0.94 | 0.68 |
143+
| `cme` | 0.0 | 0.80 | 0.69 | 0.93 | 0.79 |
144+
| `cme` | 0.0 | 0.95 | 0.78 | 0.92 | 0.84 |
145+
| `cme` | 0.3 | 0.80 | 0.70 | 0.93 | 0.80 |
146+
| `cme` | 0.5 | 0.80 | 0.70 | 0.93 | 0.80 |
147+
| `eco` | 0.0 | 0.50 | 0.39 | 0.96 | 0.56 |
148+
| `eco` | 0.0 | 0.80 | 0.66 | 0.95 | 0.78 |
149+
| `eco` | 0.0 | 0.95 | 0.76 | 0.94 | 0.84 |
150+
| `eco` | 0.3 | 0.80 | 0.67 | 0.95 | 0.78 |
151+
| `eco` | 0.5 | 0.80 | 0.69 | 0.95 | 0.80 |
152+
| `mge` | 0.0 | 0.50 | 0.62 | 0.85 | 0.72 |
153+
| `mge` | 0.0 | 0.80 | 0.77 | 0.84 | 0.80 |
154+
| `mge` | 0.0 | 0.95 | 0.82 | 0.81 | 0.81 |
155+
| `mge` | 0.3 | 0.80 | 0.77 | 0.84 | 0.80 |
156+
| `mge` | 0.5 | 0.80 | 0.78 | 0.84 | 0.81 |
157+
158+
**Reading:**
159+
- **`min_score_ratio_ko` is inert** — across all four organisms, varying it
160+
0.0 → 0.3 → 0.5 changes precision/recall by ≤0.02 (mostly 0.00). It is a
161+
magic-number knob that does effectively nothing here. (Full 0.0/0.3/0.5 × g-grid
162+
in the script output; representative rows shown.)
163+
- **`min_score_ratio_g` is the real precision lever** — 0.80 → 0.95 lifts
164+
precision ~0.07–0.10 for ~0.02 recall loss. 0.50 is clearly too loose.
165+
166+
## 4. Chosen defaults and effect
167+
168+
| parameter | RAVEN / old | new default | rationale |
169+
|---|---|---|---|
170+
| `cutoff` | 1e-50 | **1e-30** | recovers the matched tail (esp. divergent genomes); still ~22 orders above the 1e-8 noise floor |
171+
| `min_score_ratio_g` | 0.8 | **0.9** | the effective precision lever; offsets the looser cut-off |
172+
| `min_score_ratio_ko` | 0.3 | 0.3 (kept) | empirically inert; retained for RAVEN parity |
173+
174+
Old default `(1e-50, 0.3, 0.8)` vs new default `(1e-30, 0.3, 0.9)`
175+
(`min_score_ratio_ko` 0.3 ≡ 0.0 here):
176+
177+
| organism | gKO prec | gKO rec | rxn rec | rxn novel |
178+
|---|---|---|---|---|
179+
| `sce` | 0.78 → 0.76 | 0.95 → 0.96 | 0.96 → 0.96 | 157 → 137 |
180+
| `cme` | 0.70 → 0.67 | 0.93 → 0.96 | 0.94 → 0.97 | 307 → 305 |
181+
| `eco` | 0.67 → 0.67 | 0.95 → 0.97 | 0.98 → 0.99 | 198 → 173 |
182+
| `mge` | 0.77 → 0.69 | **0.84 → 0.94** | **0.87 → 0.97** | 29 → 35 |
183+
184+
The divergent minimal genome gains ~10 points of recall (the case the sequence
185+
path exists for); model organisms improve slightly and `eco` emits *fewer*
186+
unannotated reactions (the tighter gene-ratio prunes spurious multi-KO genes). The
187+
small precision dip vs annotation is dominated by extra strong homologs, not
188+
weak-hit noise.
189+
190+
## 5. Whole-model cross-validation (sanity check)
191+
192+
Full reconstruction of *S. cerevisiae* two ways, at the old defaults:
193+
194+
| | annotation path (3b.4) | HMM path (3b.5) |
195+
|---|---|---|
196+
| reactions | 1355 | 1461 |
197+
| metabolites | 1501 | 1567 |
198+
| genes | 835 | 896 |
199+
200+
Reaction recall 96.3 % (1305/1355 shared, Jaccard 0.86); gene recall 96.6 %
201+
(807/835 shared, Jaccard 0.87). The annotation path also exercises the new
202+
`organism_gene_ko.tsv.xz` artefact (K14). `hmmscan` throughput ≈ 0.1 s/query
203+
against either library on 12 threads (yeast: 6021 queries in 633 s).

0 commit comments

Comments
 (0)