From c65417ed27bd4433bd8e2de311ed9ef9e9aef40c Mon Sep 17 00:00:00 2001 From: Teo Lemane Date: Wed, 30 Sep 2026 17:22:19 +0200 Subject: [PATCH 1/4] fix(profile): fix extreme coverage (wip) --- metapang/core/profile/strains.py | 60 ++++++++++++++++++++++++-------- 1 file changed, 46 insertions(+), 14 deletions(-) diff --git a/metapang/core/profile/strains.py b/metapang/core/profile/strains.py index c69abd4..5452395 100644 --- a/metapang/core/profile/strains.py +++ b/metapang/core/profile/strains.py @@ -174,20 +174,21 @@ def union(a: int, b: int) -> None: def kstar_one_se(fold_errors: list[list[float]]) -> int: - """Pick the fewest strains whose CV error is within one SE of the minimum. - - NNLS zeros spurious strains, so the CV curve drops to the true k then plateaus - (never rises). A plain argmin chases noise to k_max, the 1-SE rule stops at the elbow. + """Sequential 1-SE stop: add a strain only while it lowers held-out error by + more than one standard error, and stop at the first increment that does not. """ means = [float(np.mean(e)) for e in fold_errors] - kbest = int(np.argmin(means)) - errs = fold_errors[kbest] - se = float(np.std(errs) / math.sqrt(len(errs))) if len(errs) > 1 else 0.0 - thresh = means[kbest] + se - for k in range(len(means)): - if means[k] <= thresh: - return k - return kbest + kstar = 0 + for k in range(len(fold_errors) - 1): + cur = fold_errors[k] + if len(cur) < 2: + break + se = float(np.std(cur) / math.sqrt(len(cur))) + if means[k] - means[k + 1] > se: + kstar = k + 1 + else: + break + return kstar def kstar_paired(fold_errors: list[list[float]], min_rel_reduction: float = 0.0) -> int: @@ -416,6 +417,24 @@ def select(self, mix: Mixture, pwg: PWGraph | None = None) -> Selection: reduced, groups = collapse_columns(mix, cfg.merge_jaccard) if not np.any(reduced.y > 0) or reduced.n_strains == 0: return Selection.empty() + + # Drop extreme-coverage families before the fit + pos = reduced.y[reduced.y > 0] + med = float(np.median(pos)) + mad = float(np.median(np.abs(pos - med))) + + fence = med + 8.0 * 1.4826 * mad if mad > 0 else 20.0 * med + keep = reduced.y <= fence + if not bool(keep.all()): + fl = reduced.family_labels + reduced = Mixture( + m=reduced.m[keep], + y=reduced.y[keep], + strain_labels=reduced.strain_labels, + family_labels=[f for f, k in zip(fl, keep, strict=True) if k] + if fl + else None, + ) return self._select_cv(reduced, groups, pwg) @@ -588,12 +607,25 @@ def refine_abundances(selection: Selection, components: list[GeneSet]) -> Select M[i, cid] = 1.0 x, _ = nnls(M, y) + + # Drop strains the refit zeroed out + keep = [i for i in range(n) if x[i] > 1e-9] + if not keep: + return Selection.empty() + if len(keep) < n: + x, _ = nnls(M[:, keep], y) + indices = [selection.indices[i] for i in keep] + labels = [selection.labels[i] for i in keep] + else: + indices = selection.indices + labels = selection.labels + total = float(np.sum(x)) ra = x / total if total > 0 else x return Selection( - indices=selection.indices, - labels=selection.labels, + indices=indices, + labels=labels, x=x, ra=ra, residual=selection.residual, From 213dee4bfcddedf50ac9738290928e5d290689c2 Mon Sep 17 00:00:00 2001 From: Teo Lemane Date: Thu, 1 Oct 2026 09:29:45 +0200 Subject: [PATCH 2/4] fix(profile): use poisson deviance (wip) --- metapang/core/profile/strains.py | 52 ++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 12 deletions(-) diff --git a/metapang/core/profile/strains.py b/metapang/core/profile/strains.py index 5452395..f904d6a 100644 --- a/metapang/core/profile/strains.py +++ b/metapang/core/profile/strains.py @@ -249,6 +249,16 @@ def _nnls_via_gram( return x, max(rss, 0.0) +def _poisson_deviance(y: NDArray[np.float64], mu: NDArray[np.float64]) -> float: + """Poisson deviance""" + y = np.asarray(y, dtype=float) + mu = np.asarray(mu, dtype=float) + term = np.zeros_like(y) + pos = y > 0 + term[pos] = y[pos] * np.log(y[pos] / mu[pos]) + return float(2.0 * np.sum(term - (y - mu))) + + class GreedyStrainSelector: """Selects strains by greedy forward NNLS with cross-validated stopping.""" @@ -266,6 +276,10 @@ def __init__( self._fold: list[ tuple[NDArray[np.float64], NDArray[np.float64], float, int] ] = [] + self._M: NDArray[np.float64] = np.zeros((0, 0)) + self._y: NDArray[np.float64] = np.zeros(0) + self._fold_id: NDArray[np.int_] = np.zeros(0, dtype=int) + self._b0: float = 1e-6 def _solve_subset(self, M, y, cols, labels): return _nnls_via_gram(self._G, self._b, self._yty, cols) @@ -284,20 +298,27 @@ def _pick_best(self, M, y, selected, remaining, labels): return tied[0] def _cv_fold_errors(self, cols, n_folds: int) -> list[float]: - """Per-fold cross-validated MSE for the model using `cols`: fit on training folds, predict the held-out fold.""" + """Per-fold cross-validated Poisson deviance for the model using `cols`: fit + on the training folds, then score the held-out fold with + the count-appropriate loss. + """ cols = list(cols) errs: list[float] = [] for f in range(n_folds): Gf, bf, ytf, cnt = self._fold[f] if cnt == 0 or cnt == self._n_fam: continue + rows = self._fold_id == f + yf = self._y[rows] if not cols: - errs.append(ytf / cnt) - continue - xx, _ = _nnls_via_gram(self._G - Gf, self._b - bf, self._yty - ytf, cols) - Gfc = Gf[np.ix_(cols, cols)] - test_rss = float(xx @ Gfc @ xx - 2.0 * xx @ bf[cols] + ytf) - errs.append(max(test_rss, 0.0) / cnt) + mu = np.zeros(yf.shape[0]) + else: + xx, _ = _nnls_via_gram( + self._G - Gf, self._b - bf, self._yty - ytf, cols + ) + mu = self._M[np.ix_(rows, cols)] @ xx + dev = _poisson_deviance(yf, mu + self._b0) + errs.append(dev / cnt) return errs def _finalize(self, reduced, sel_cols, groups, trace) -> Selection: @@ -375,9 +396,17 @@ def _select_cv(self, reduced, groups, pwg) -> Selection: self._yty = float(y @ y) self._n_fam = n_fam + # kept so the CV scorer can materialise per-family held-out predictions for + # the Poisson deviance + self._M = M + self._y = y + pos = y[y > 0] + self._b0 = 1e-2 * float(np.median(pos)) if pos.size else 1e-6 + # Per-fold Gram pieces (same fixed partition/seed as before) so CV folds # reuse the precomputed Gram: G_train = G - G_fold, etc. fold = np.random.default_rng(0).integers(0, cfg.cv_folds, size=n_fam) + self._fold_id = fold self._fold = [] for f in range(cfg.cv_folds): rows = fold == f @@ -418,12 +447,11 @@ def select(self, mix: Mixture, pwg: PWGraph | None = None) -> Selection: if not np.any(reduced.y > 0) or reduced.n_strains == 0: return Selection.empty() - # Drop extreme-coverage families before the fit + # Drop extreme-coverage families before the fit, using a Tukey fence pos = reduced.y[reduced.y > 0] - med = float(np.median(pos)) - mad = float(np.median(np.abs(pos - med))) - - fence = med + 8.0 * 1.4826 * mad if mad > 0 else 20.0 * med + q1, med, q3 = (float(v) for v in np.percentile(pos, [25, 50, 75])) + iqr = q3 - q1 + fence = q3 + 3.0 * iqr if iqr > 0 else 20.0 * med keep = reduced.y <= fence if not bool(keep.all()): fl = reduced.family_labels From 3e4f33e1c23273f286f6b4f864ad45e44afed482 Mon Sep 17 00:00:00 2001 From: Teo Lemane Date: Thu, 1 Oct 2026 09:32:27 +0200 Subject: [PATCH 3/4] feat(profile): output species detection scores --- metapang/app/profile/__init__.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/metapang/app/profile/__init__.py b/metapang/app/profile/__init__.py index 9acc008..4f638dd 100644 --- a/metapang/app/profile/__init__.py +++ b/metapang/app/profile/__init__.py @@ -255,6 +255,18 @@ def gather_candidates(): candidates.append(species) mp_log.trace(f" candidate species '{species}' (from genome '{name}')") + detection_file = output_directory / f"{sample_name}.detection.tsv" + with open(detection_file, "w") as det: + det.write("genome\tspecies\tf_query\tf_match\tcum_f_query\n") + cum = 0.0 + for name, f_query, f_match in gather_results: + cum += float(f_query) + det.write( + f"{name}\t{name.split('@')[0]}\t" + f"{float(f_query):.6f}\t{float(f_match):.6f}\t{cum:.6f}\n" + ) + mp_log.info(f"Wrote species detection scores to '{detection_file.name}'") + if candidates: mp_log.info(f"{len(candidates)} candidate species: {', '.join(candidates)}") else: From 3736ee2f74403eaa85085a4d8b00f73300ce9f52 Mon Sep 17 00:00:00 2001 From: Teo Lemane Date: Wed, 21 Oct 2026 11:19:15 +0200 Subject: [PATCH 4/4] fix(core): batch signature threading fix(core): batch signature threading --- metapang/core/index/index.py | 197 +++++++++++++++++++++++------------ 1 file changed, 129 insertions(+), 68 deletions(-) diff --git a/metapang/core/index/index.py b/metapang/core/index/index.py index 8dbd192..bdb0278 100644 --- a/metapang/core/index/index.py +++ b/metapang/core/index/index.py @@ -1,6 +1,7 @@ import json import os import signal +import threading from collections.abc import Iterable from concurrent.futures import ProcessPoolExecutor, as_completed from dataclasses import asdict, dataclass, field @@ -246,27 +247,38 @@ def pangenome_construct( @staticmethod def file_batch_signature( - file: "Path | list[Path]", - k: int, - scaled: int, - n: int, - force: bool, - nb_cores: int, - queue_size: int, - batch_size: int, - with_log: bool = True, - ) -> SourmashSignature: - """Compute one Sourmash signature for one or several files, in parallel batches. - - Several files (e.g. paired-end or split reads) are folded into a single MinHash, - so the result is one signature for the whole sample. - """ + file, k, scaled, n, force, nb_cores, queue_size, batch_size, with_log=True + ): files = [file] if isinstance(file, (str, Path)) else list(file) - man = MinHashManager(ctx=_MP_CONTEXT) - man.register("MinHash", MinHash) - man.start() - m = man.MinHash(n=n, ksize=k, scaled=scaled) # type: ignore[attr-defined] + result = MinHash(n=n, ksize=k, scaled=scaled) + merge_lock = threading.Lock() # callbacks run in a thread of the main process + errors: list[BaseException] = [] + + def _merge(fut): + try: + mh = fut.result() + except BaseException as e: + errors.append(e) + return + with merge_lock: + result.merge(mh) + + def _iter_sequences(): + for f in files: + if is_fastx(f): + for record in Fastx(f): + yield record[1] + else: + with zopen(f, "rb") as file_io: + for record in gb_io.iter(file_io): + yield record.sequence.decode("utf-8") + + def _submit(executor, batch): + fut = executor.submit( + IndexBuilder.sequence_batch_minhash, batch, k, n, scaled, force + ) + fut.add_done_callback(_merge) with BoundedProcessPoolExecutor( queue_size, @@ -275,60 +287,109 @@ def file_batch_signature( initializer=_ignore_sigint, ) as executor: try: - sequences = [] - i = 0 - for f in files: - if is_fastx(f): - for record in Fastx(f): - sequences.append(record[1]) - i += 1 - if i % batch_size == 0: - executor.submit( - IndexBuilder.add_sequences_minhash, - m, - sequences, - k, - n, - scaled, - force, - ) - sequences = [] - else: - with zopen(f, "rb") as file_io: - for record in gb_io.iter(file_io): - sequences.append(record.sequence.decode("utf-8")) - i += 1 - if i % batch_size == 0: - executor.submit( - IndexBuilder.add_sequences_minhash, - m, - sequences, - k, - n, - scaled, - force, - ) - sequences = [] - - if sequences: - executor.submit( - IndexBuilder.add_sequences_minhash, - m, - sequences, - k, - n, - scaled, - force, - ) - + batch = [] + for seq in _iter_sequences(): + batch.append(seq) + if len(batch) == batch_size: + _submit(executor, batch) + batch = [] + if batch: + _submit(executor, batch) except KeyboardInterrupt: mp_log_error(with_log, "Interrupted by user (Ctrl+C)") executor.shutdown(wait=False) os._exit(1) - signature = SourmashSignature(m.copy(), name=files[0].stem) - man.shutdown() - return signature + if errors: + raise errors[0] + + return SourmashSignature(result, name=Path(files[0]).stem) + + # @staticmethod + # def file_batch_signature( + # file: "Path | list[Path]", + # k: int, + # scaled: int, + # n: int, + # force: bool, + # nb_cores: int, + # queue_size: int, + # batch_size: int, + # with_log: bool = True, + # ) -> SourmashSignature: + # """Compute one Sourmash signature for one or several files, in parallel batches. + + # Several files (e.g. paired-end or split reads) are folded into a single MinHash, + # so the result is one signature for the whole sample. + # """ + # files = [file] if isinstance(file, (str, Path)) else list(file) + + # man = MinHashManager(ctx=_MP_CONTEXT) + # man.register("MinHash", MinHash) + # man.start() + # m = man.MinHash(n=n, ksize=k, scaled=scaled) # type: ignore[attr-defined] + + # with BoundedProcessPoolExecutor( + # queue_size, + # max_workers=nb_cores, + # mp_context=_MP_CONTEXT, + # initializer=_ignore_sigint, + # ) as executor: + # try: + # sequences = [] + # i = 0 + # for f in files: + # if is_fastx(f): + # for record in Fastx(f): + # sequences.append(record[1]) + # i += 1 + # if i % batch_size == 0: + # executor.submit( + # IndexBuilder.add_sequences_minhash, + # m, + # sequences, + # k, + # n, + # scaled, + # force, + # ) + # sequences = [] + # else: + # with zopen(f, "rb") as file_io: + # for record in gb_io.iter(file_io): + # sequences.append(record.sequence.decode("utf-8")) + # i += 1 + # if i % batch_size == 0: + # executor.submit( + # IndexBuilder.add_sequences_minhash, + # m, + # sequences, + # k, + # n, + # scaled, + # force, + # ) + # sequences = [] + + # if sequences: + # executor.submit( + # IndexBuilder.add_sequences_minhash, + # m, + # sequences, + # k, + # n, + # scaled, + # force, + # ) + + # except KeyboardInterrupt: + # mp_log_error(with_log, "Interrupted by user (Ctrl+C)") + # executor.shutdown(wait=False) + # os._exit(1) + + # signature = SourmashSignature(m.copy(), name=files[0].stem) + # man.shutdown() + # return signature def construct(self, nb_cores: int, with_log: bool = True) -> IndexInfo: """Construct and save the index, returning its IndexInfo."""