Skip to content

Commit 97d3339

Browse files
committed
Merge dev: gz inspect.txt support + unmapped-validation-taxid warning
2 parents 1fc797a + 6614ae3 commit 97d3339

4 files changed

Lines changed: 150 additions & 12 deletions

File tree

nanometa_live/app/tabs/kraken2_helpers.py

Lines changed: 30 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -203,25 +203,41 @@ def load_kraken2_taxonomy(kraken_db_path: str) -> dict:
203203
Dict mapping taxid -> parent_taxid. Returns empty dict if the
204204
inspect.txt file is missing or cannot be read.
205205
"""
206+
import gzip
206207
import os
207208

208209
if not kraken_db_path:
209210
return {}
210211

211-
inspect_path = os.path.join(kraken_db_path, "inspect.txt")
212-
if inspect_path in _TAXONOMY_CACHE:
213-
return _TAXONOMY_CACHE[inspect_path]
214-
215-
if not os.path.exists(inspect_path):
216-
logging.debug(f"Kraken2 inspect.txt not found at {inspect_path}")
217-
_TAXONOMY_CACHE[inspect_path] = {}
212+
# Cache on the database directory so a gz-only DB is cached too (the cache
213+
# key used to be the plain inspect.txt path).
214+
if kraken_db_path in _TAXONOMY_CACHE:
215+
return _TAXONOMY_CACHE[kraken_db_path]
216+
217+
# Prefer a plain inspect.txt; fall back to a gzipped inspect.txt.gz. Some
218+
# Kraken2 builds (e.g. GTDB-derived / size-conscious DBs) ship only the
219+
# compressed form. The DB indexer already reads the .gz
220+
# (database_indexer.build_from_inspect_gz); mirror that here so the
221+
# authoritative-taxonomy correction for Sankey/Sunburst is not silently
222+
# disabled on those databases.
223+
plain_path = os.path.join(kraken_db_path, "inspect.txt")
224+
gz_path = os.path.join(kraken_db_path, "inspect.txt.gz")
225+
if os.path.exists(plain_path):
226+
inspect_path = plain_path
227+
_opener = lambda p: open(p)
228+
elif os.path.exists(gz_path):
229+
inspect_path = gz_path
230+
_opener = lambda p: gzip.open(p, "rt", encoding="utf-8")
231+
else:
232+
logging.debug(f"Kraken2 inspect.txt[.gz] not found in {kraken_db_path}")
233+
_TAXONOMY_CACHE[kraken_db_path] = {}
218234
return {}
219235

220236
taxid_to_parent: dict = {}
221237
indent_stack = [] # list of (indent, taxid)
222238

223239
try:
224-
with open(inspect_path) as f:
240+
with _opener(inspect_path) as f:
225241
for line in f:
226242
if line.startswith("#"):
227243
continue
@@ -241,16 +257,18 @@ def load_kraken2_taxonomy(kraken_db_path: str) -> dict:
241257
parent = indent_stack[-1][1] if indent_stack else 0
242258
taxid_to_parent[taxid] = parent
243259
indent_stack.append((indent, taxid))
244-
except OSError as exc:
245-
logging.warning(f"Failed to read Kraken2 inspect.txt: {exc}")
246-
_TAXONOMY_CACHE[inspect_path] = {}
260+
except (OSError, EOFError) as exc:
261+
# EOFError / gzip.BadGzipFile (an OSError subclass) cover a truncated or
262+
# malformed .gz; a plain-text OSError covers an unreadable inspect.txt.
263+
logging.warning(f"Failed to read Kraken2 {os.path.basename(inspect_path)}: {exc}")
264+
_TAXONOMY_CACHE[kraken_db_path] = {}
247265
return {}
248266

249267
logging.info(
250268
f"Loaded Kraken2 taxonomy from {inspect_path}: "
251269
f"{len(taxid_to_parent)} taxa"
252270
)
253-
_TAXONOMY_CACHE[inspect_path] = taxid_to_parent
271+
_TAXONOMY_CACHE[kraken_db_path] = taxid_to_parent
254272
return taxid_to_parent
255273

256274

nanometa_live/core/config/parameter_mapping.py

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,12 @@ def get_validation_species_from_watchlist(
148148

149149
species_list = []
150150
genome_paths = []
151+
# Track entries whose Kraken2 taxid fell back to the raw NCBI taxid
152+
# because nothing mapped them. On a GTDB/custom-taxid database (where the
153+
# NCBI taxid is absent), read extraction by that taxid yields ZERO reads
154+
# and validation silently produces nothing -- the failure mode found
155+
# auditing a custom-taxid DB. We surface it as a warning below.
156+
unmapped_names = []
151157

152158
for entry in entries:
153159
ncbi_taxid = getattr(entry, 'taxid', 0)
@@ -159,9 +165,11 @@ def get_validation_species_from_watchlist(
159165
# so the operator does not have to run "Scan Database"; otherwise
160166
# fall back to the auto-mapping collection, then to the NCBI taxid.
161167
kraken_taxid = ncbi_taxid # Default to NCBI taxid
168+
mapped = False
162169
explicit_db_taxid = getattr(entry, 'db_taxid', None)
163170
if explicit_db_taxid:
164171
kraken_taxid = explicit_db_taxid
172+
mapped = True
165173
logging.debug(
166174
f"Using explicit db_taxid {explicit_db_taxid} for "
167175
f"{getattr(entry, 'name', '')}"
@@ -170,10 +178,14 @@ def get_validation_species_from_watchlist(
170178
db_taxid = mapping_collection.get_db_taxid(ncbi_taxid)
171179
if db_taxid:
172180
kraken_taxid = db_taxid
181+
mapped = True
173182
logging.debug(
174183
f"Mapped NCBI {ncbi_taxid} -> Kraken2 {db_taxid} for {getattr(entry, 'name', '')}"
175184
)
176185

186+
if not mapped:
187+
unmapped_names.append(getattr(entry, 'name', '') or str(ncbi_taxid))
188+
177189
species_info = {
178190
'taxid': ncbi_taxid,
179191
'kraken_taxid': kraken_taxid,
@@ -190,6 +202,31 @@ def get_validation_species_from_watchlist(
190202
f"Found {len(species_list)} enabled watchlist species, "
191203
f"{len(genome_paths)} with downloaded genomes"
192204
)
205+
206+
# Warn loudly when validation taxids fell back to raw NCBI taxids: on a
207+
# custom/GTDB database this silently extracts 0 reads. Distinguish "never
208+
# scanned" (no mapping at all) from "scanned but these did not map".
209+
if unmapped_names:
210+
preview = ", ".join(unmapped_names[:5]) + (
211+
f" (+{len(unmapped_names) - 5} more)" if len(unmapped_names) > 5 else "")
212+
if mapping_collection is None:
213+
logging.warning(
214+
"No Kraken2 taxid mapping has been generated for this database, "
215+
"so validation will use raw NCBI taxids (%d species: %s). If this "
216+
"is a GTDB or custom database whose taxids differ from NCBI, read "
217+
"extraction will find ZERO reads and validation will produce "
218+
"nothing. Run 'Scan Database' / 'Verify Taxonomy IDs' in the "
219+
"Watchlist & Preparation tab before starting.",
220+
len(unmapped_names), preview,
221+
)
222+
else:
223+
logging.warning(
224+
"%d watchlist species could not be mapped to this Kraken2 "
225+
"database (%s); validation will fall back to their NCBI taxids "
226+
"and may find no reads if the database does not contain them.",
227+
len(unmapped_names), preview,
228+
)
229+
193230
return species_list, genome_paths
194231

195232
except (ImportError, AttributeError) as e:

tests/test_kraken2_helpers.py

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,30 @@ def test_missing_inspect_returns_empty(self, tmp_path):
6363
def test_empty_path_returns_empty(self):
6464
assert load_kraken2_taxonomy("") == {}
6565

66+
def test_parses_gzipped_inspect(self, tmp_path):
67+
# GTDB-derived / size-conscious DBs ship only inspect.txt.gz. The
68+
# authoritative-taxonomy correction must still load (was silently empty).
69+
import gzip
70+
content = (
71+
"100.00\t1000\t0\tR\t1\troot\n"
72+
"90.00\t900\t0\tD\t2\t Bacteria\n"
73+
"50.00\t500\t0\tG\t561\t Escherichia\n"
74+
"25.00\t250\t250\tS\t562\t Escherichia coli\n"
75+
)
76+
with gzip.open(tmp_path / "inspect.txt.gz", "wt", encoding="utf-8") as fh:
77+
fh.write(content)
78+
mapping = load_kraken2_taxonomy(str(tmp_path))
79+
assert mapping == {1: 0, 2: 1, 561: 2, 562: 561}
80+
81+
def test_plain_inspect_preferred_over_gz(self, tmp_path):
82+
import gzip
83+
(tmp_path / "inspect.txt").write_text(
84+
"100.00\t1000\t0\tR\t1\troot\n90.00\t900\t0\tD\t2\t Bacteria\n"
85+
)
86+
with gzip.open(tmp_path / "inspect.txt.gz", "wt", encoding="utf-8") as fh:
87+
fh.write("bogus\tnot\tparsed\n") # would corrupt if gz were read
88+
assert load_kraken2_taxonomy(str(tmp_path)) == {1: 0, 2: 1}
89+
6690

6791
class TestApplyAuthoritativeTaxonomy:
6892
def test_replaces_parent_taxid_from_mapping(self):

tests/test_parameter_mapping_extra.py

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -149,6 +149,65 @@ def test_empty_watchlist(self):
149149
assert get_validation_species({}) == ([], [])
150150

151151

152+
class TestUnmappedValidationTaxidWarning:
153+
"""A custom/GTDB DB whose taxids differ from NCBI silently extracts 0 reads
154+
when the watchlist was never mapped. get_validation_species_from_watchlist
155+
must WARN at that fallback instead of failing silently."""
156+
157+
def _entry(self, taxid, name, db_taxid=None):
158+
e = MagicMock()
159+
e.taxid = taxid
160+
e.name = name
161+
e.db_taxid = db_taxid
162+
e.names_alt = []
163+
return e
164+
165+
def _run(self, entries, mapping_collection):
166+
wm = MagicMock()
167+
wm._loaded = True
168+
wm.get_active_entries.return_value = {e.taxid: e for e in entries}
169+
gm = MagicMock()
170+
gm.get_genome_path.return_value = None
171+
with patch.object(pm, "get_watchlist_manager", return_value=wm), \
172+
patch.object(pm, "get_genome_manager", return_value=gm), \
173+
patch("nanometa_live.core.taxonomy.taxid_mapping.get_mapping_collection",
174+
return_value=mapping_collection):
175+
return pm.get_validation_species_from_watchlist({"kraken_db": "/db"})
176+
177+
def test_no_mapping_warns_run_scan(self, caplog):
178+
import logging
179+
with caplog.at_level(logging.WARNING):
180+
species, _ = self._run([self._entry(263, "Francisella tularensis")], None)
181+
assert species and species[0]["kraken_taxid"] == 263 # raw NCBI fallback
182+
assert any("Scan Database" in r.message for r in caplog.records)
183+
184+
def test_scanned_all_mapped_no_warning(self, caplog):
185+
import logging
186+
coll = MagicMock()
187+
coll.get_db_taxid.return_value = 4007169 # 263 -> custom DB taxid
188+
with caplog.at_level(logging.WARNING):
189+
species, _ = self._run([self._entry(263, "Francisella tularensis")], coll)
190+
assert species[0]["kraken_taxid"] == 4007169
191+
assert not any("could not be mapped" in r.message or "Scan Database" in r.message
192+
for r in caplog.records)
193+
194+
def test_scanned_partial_warns_specific(self, caplog):
195+
import logging
196+
coll = MagicMock()
197+
coll.get_db_taxid.side_effect = lambda t: 4007169 if t == 263 else None
198+
with caplog.at_level(logging.WARNING):
199+
self._run([self._entry(263, "F. tularensis"),
200+
self._entry(1392, "B. anthracis")], coll)
201+
assert any("could not be mapped" in r.message for r in caplog.records)
202+
203+
def test_explicit_db_taxid_not_warned(self, caplog):
204+
import logging
205+
with caplog.at_level(logging.WARNING):
206+
species, _ = self._run([self._entry(263, "F. tularensis", db_taxid=4007169)], None)
207+
assert species[0]["kraken_taxid"] == 4007169
208+
assert not any("Scan Database" in r.message for r in caplog.records)
209+
210+
152211
class TestPathogenGenomesLocation:
153212
"""Operator feedback #2: archive/rerun crashed with 'No such file:
154213
.../validation/pathogen_genomes.json'. The launch input must live OUTSIDE

0 commit comments

Comments
 (0)