From 1fe56bf334c715a4c78eb738c3b22adde2888843 Mon Sep 17 00:00:00 2001 From: marco cominelli Date: Thu, 26 Mar 2026 12:17:41 +0000 Subject: [PATCH 1/6] Add filter for unreliable features and additional convergence criterion --- irescue/count.py | 43 ++++++++++++++++++++++++++++++++----------- irescue/em.py | 20 +++++++++++--------- irescue/main.py | 21 +++++++++++++++++++++ 3 files changed, 64 insertions(+), 20 deletions(-) diff --git a/irescue/count.py b/irescue/count.py index 2f544f2..08de4da 100644 --- a/irescue/count.py +++ b/irescue/count.py @@ -94,7 +94,7 @@ def parse_maps(maps_file, feature_index): def compute_cell_counts( - equivalence_classes, features_index, max_iters, tolerance, dumpEC, no_umi + equivalence_classes, features_index, max_iters, tolerance, dumpEC, no_umi, exclude_unreliable_features, convergence_criterion ): """ Calculate TE counts of a single cell, given a list of equivalence classes. @@ -259,18 +259,35 @@ def compute_cell_counts( em_array = em_array.tocsr() # save an array with features > 0, as in em_array order - tokeep = np.flatnonzero(em_array.sum(axis=0)) + if not exclude_unreliable_features: + tokeep = np.flatnonzero(em_array.sum(axis=0)) + else: + # keep features supported by at least 2 multimapping + tokeep1 = np.where((em_array.sum(axis=0) >= 2).A1)[0] + # rescue features supported by only 1 multimapping but also by (at least) 1 uniquely mapping + tokeep2 = np.intersect1d( + np.array(list(counts.keys())), + np.where((em_array.sum(axis=0) == 1).A1)[0] + ) + tokeep = np.union1d(tokeep1, tokeep2) + # remove unmapped features from em_array em_array = em_array[:, tokeep] - # run EM - em_counts, em_stats = run_em( - em_array, cycles=max_iters, tolerance=tolerance - ) - em_counts = em_counts * em_array.shape[0] - for i, c in zip(tokeep + 1, em_counts): - if c > 0: - counts[i] += c + # removing some features may yield empty rows (not necessary if exclude_unreliable_features is disabled) + if exclude_unreliable_features: + em_array = em_array[(em_array.sum(axis=1)>0).A1, :] + + if em_array.shape[1] > 0: + # run EM + em_counts, em_stats = run_em( + em_array, cycles=max_iters, tolerance=tolerance, convergence_criterion=convergence_criterion + ) + em_counts = em_counts * em_array.shape[0] + + for i, c in zip(tokeep + 1, em_counts): + if c > 0: + counts[i] += c return dict(counts), dump, em_stats @@ -292,6 +309,8 @@ def run_count( features_index, tmpdir, no_umi, + exclude_unreliable_features, + convergence_criterion, dumpEC, max_iters, tolerance, @@ -323,11 +342,13 @@ def run_count( tolerance=tolerance, dumpEC=dumpEC, no_umi=no_umi, + exclude_unreliable_features=exclude_unreliable_features, + convergence_criterion=convergence_criterion ) writerr( f"[{taskn}] Write cell {cellidx} ({cellbarcode.decode()}). " f"EM cycles: {em_stats[0]}. Converged: {em_stats[1]}. " - f"Log likelihood: {em_stats[2]}. Increment: {em_stats[3]}.", + f"Log likelihood: {em_stats[2] if convergence_criterion=='likelihood' else 'Not computed because of convergence criterion choice'}. Increment: {em_stats[3]}.", level=1, send=verbose, ) diff --git a/irescue/em.py b/irescue/em.py index 32c94e7..6a48b85 100644 --- a/irescue/em.py +++ b/irescue/em.py @@ -29,7 +29,7 @@ def log_likelihood(matrix, counts): return log_likelihood -def run_em(matrix, cycles=100, tolerance=1e-4): +def run_em(matrix, cycles=100, tolerance=1e-4, convergence_criterion="likelihood"): """ Run Expectation-Maximization (EM) algorithm to redistribute read counts across a set of features. @@ -42,6 +42,8 @@ def run_em(matrix, cycles=100, tolerance=1e-4): Number of EM cycles. tolerance : float Tolerance threshold of log-likelihood difference to infer convergence. + convergence_criterion : str + Criterion to determine convergence: "likelihood" or "parameters". Returns ------- @@ -59,8 +61,8 @@ def run_em(matrix, cycles=100, tolerance=1e-4): nFeatures = matrix.shape[1] counts = np.full(shape=nFeatures, fill_value=1 / nFeatures) - # Initial log-likelihood - prev_loglik = log_likelihood(matrix, counts) + # Initial log-likelihood (or initial counts) + prev = log_likelihood(matrix, counts) if convergence_criterion=="likelihood" else counts converged = False curr_cycle = 0 @@ -71,15 +73,15 @@ def run_em(matrix, cycles=100, tolerance=1e-4): e_matrix = e_step(matrix=matrix, counts=counts) counts = m_step(matrix=e_matrix) - # Compute the new log-likelihood - loglik = log_likelihood(matrix, counts) + # Compute the new log-likelihood (depending on convergence criterion) + curr = log_likelihood(matrix, counts) if convergence_criterion=="likelihood" else counts # Check for convergence - loglikdiff = loglik - prev_loglik - if np.abs(loglikdiff) < tolerance: + diff = np.abs(curr-prev) if convergence_criterion=="likelihood" else np.abs(curr-prev).sum() + if diff < tolerance: converged = True break - prev_loglik = loglik + prev = curr - return counts, (curr_cycle, converged, loglik, loglikdiff) + return counts, (curr_cycle, converged, curr, diff) diff --git a/irescue/main.py b/irescue/main.py index 226197e..24b46a5 100644 --- a/irescue/main.py +++ b/irescue/main.py @@ -181,6 +181,25 @@ def parseArguments(): "(Default: %(default)s)." ), ) + parser.add_argument( + "--exclude-unreliable-features", + action="store_true", + help=( + "Exclude features supported by only 1 multimapping read " + "(Default: %(default)s)." + ), + ) + parser.add_argument( + "--convergence-criterion", + type=str, + metavar="STR", + choices=["likelihood", "parameters"], + default="likelihood", + help=( + "Criterion to define convergence. " + "One of: likelihood, parameters (Default: %(default)s)." + ), + ) parser.add_argument( "--dump-ec", action="store_true", @@ -379,6 +398,8 @@ def main(): feature_index, dirs["tmp"], args.no_umi, + args.exclude_unreliable_features, + args.convergence_criterion, args.dump_ec, args.max_iters, args.tolerance, From b2ea8b9675fd13386659e9dcd1228fb24eb1c2bf Mon Sep 17 00:00:00 2001 From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com> Date: Wed, 15 Apr 2026 14:42:32 +0200 Subject: [PATCH 2/6] fix base index mismatch between `em_array` and `counts` --- irescue/count.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/irescue/count.py b/irescue/count.py index 08de4da..2260fd2 100644 --- a/irescue/count.py +++ b/irescue/count.py @@ -262,11 +262,13 @@ def compute_cell_counts( if not exclude_unreliable_features: tokeep = np.flatnonzero(em_array.sum(axis=0)) else: - # keep features supported by at least 2 multimapping + # only keep features supported by >=2 multimapping reads, + # or >=1 multimapping and >=1 uniquely mapping tokeep1 = np.where((em_array.sum(axis=0) >= 2).A1)[0] - # rescue features supported by only 1 multimapping but also by (at least) 1 uniquely mapping tokeep2 = np.intersect1d( - np.array(list(counts.keys())), + # (-1 because of 0-based indexing of em_array + # against 1-based indexing of features) + np.array(list(counts.keys())) - 1, np.where((em_array.sum(axis=0) == 1).A1)[0] ) tokeep = np.union1d(tokeep1, tokeep2) @@ -274,7 +276,8 @@ def compute_cell_counts( # remove unmapped features from em_array em_array = em_array[:, tokeep] - # removing some features may yield empty rows (not necessary if exclude_unreliable_features is disabled) + # removing some features may yield empty rows + # (not necessary if exclude_low_support is disabled) if exclude_unreliable_features: em_array = em_array[(em_array.sum(axis=1)>0).A1, :] @@ -285,6 +288,10 @@ def compute_cell_counts( ) em_counts = em_counts * em_array.shape[0] + # add EM-optimized counts to uniquely mapped counts + # (add +1 to features to multimapped features to keep + # because of 0-based indexing of em_array against + # 1-based indexing of features) for i, c in zip(tokeep + 1, em_counts): if c > 0: counts[i] += c From 076149a2823e74b3671871cc99337c241c7168d1 Mon Sep 17 00:00:00 2001 From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com> Date: Wed, 15 Apr 2026 14:44:22 +0200 Subject: [PATCH 3/6] rename --exclude-unreliable-features to --exclude-low-support --- irescue/count.py | 10 +++++----- irescue/main.py | 23 ++++++++++++----------- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/irescue/count.py b/irescue/count.py index 2260fd2..94274dc 100644 --- a/irescue/count.py +++ b/irescue/count.py @@ -94,7 +94,7 @@ def parse_maps(maps_file, feature_index): def compute_cell_counts( - equivalence_classes, features_index, max_iters, tolerance, dumpEC, no_umi, exclude_unreliable_features, convergence_criterion + equivalence_classes, features_index, max_iters, tolerance, dumpEC, no_umi, exclude_low_support, convergence_criterion ): """ Calculate TE counts of a single cell, given a list of equivalence classes. @@ -259,7 +259,7 @@ def compute_cell_counts( em_array = em_array.tocsr() # save an array with features > 0, as in em_array order - if not exclude_unreliable_features: + if not exclude_low_support: tokeep = np.flatnonzero(em_array.sum(axis=0)) else: # only keep features supported by >=2 multimapping reads, @@ -278,7 +278,7 @@ def compute_cell_counts( # removing some features may yield empty rows # (not necessary if exclude_low_support is disabled) - if exclude_unreliable_features: + if exclude_low_support: em_array = em_array[(em_array.sum(axis=1)>0).A1, :] if em_array.shape[1] > 0: @@ -316,7 +316,7 @@ def run_count( features_index, tmpdir, no_umi, - exclude_unreliable_features, + exclude_low_support, convergence_criterion, dumpEC, max_iters, @@ -349,7 +349,7 @@ def run_count( tolerance=tolerance, dumpEC=dumpEC, no_umi=no_umi, - exclude_unreliable_features=exclude_unreliable_features, + exclude_low_support=exclude_low_support, convergence_criterion=convergence_criterion ) writerr( diff --git a/irescue/main.py b/irescue/main.py index 24b46a5..8bf3bb1 100644 --- a/irescue/main.py +++ b/irescue/main.py @@ -177,15 +177,8 @@ def parseArguments(): metavar="FLOAT", default=1e-4, help=( - "Log-likelihood change below which convergence is assumed " - "(Default: %(default)s)." - ), - ) - parser.add_argument( - "--exclude-unreliable-features", - action="store_true", - help=( - "Exclude features supported by only 1 multimapping read " + "Change between EM iterations below which convergence is assumed, " + "calculated on --convergence-criterion." "(Default: %(default)s)." ), ) @@ -197,7 +190,15 @@ def parseArguments(): default="likelihood", help=( "Criterion to define convergence. " - "One of: likelihood, parameters (Default: %(default)s)." + "One of: %(choices)s. (Default: %(default)s)." + ), + ) + parser.add_argument( + "--exclude-low-support", + action="store_true", + help=( + "Exclude features supported by only 1 multimapping read. " + "This should improve performance without significantly affecting results." ), ) parser.add_argument( @@ -398,7 +399,7 @@ def main(): feature_index, dirs["tmp"], args.no_umi, - args.exclude_unreliable_features, + args.exclude_low_support, args.convergence_criterion, args.dump_ec, args.max_iters, From 9d71acc6950b825a928aa6d7c2ca9f36b6dad6ae Mon Sep 17 00:00:00 2001 From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com> Date: Wed, 15 Apr 2026 14:45:36 +0200 Subject: [PATCH 4/6] better choices documentation for --converge-criterion --- irescue/em.py | 19 +++++++++++++++---- irescue/main.py | 2 +- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/irescue/em.py b/irescue/em.py index 6a48b85..9b0187a 100644 --- a/irescue/em.py +++ b/irescue/em.py @@ -43,7 +43,9 @@ def run_em(matrix, cycles=100, tolerance=1e-4, convergence_criterion="likelihood tolerance : float Tolerance threshold of log-likelihood difference to infer convergence. convergence_criterion : str - Criterion to determine convergence: "likelihood" or "parameters". + Criterion to determine convergence: + "likelihood": log-likelihood change < tolerance. + "estimates": feature abundances change < tolerance. Returns ------- @@ -62,7 +64,10 @@ def run_em(matrix, cycles=100, tolerance=1e-4, convergence_criterion="likelihood counts = np.full(shape=nFeatures, fill_value=1 / nFeatures) # Initial log-likelihood (or initial counts) - prev = log_likelihood(matrix, counts) if convergence_criterion=="likelihood" else counts + if convergence_criterion=="likelihood": + prev = log_likelihood(matrix, counts) + else: + prev = counts converged = False curr_cycle = 0 @@ -74,10 +79,16 @@ def run_em(matrix, cycles=100, tolerance=1e-4, convergence_criterion="likelihood counts = m_step(matrix=e_matrix) # Compute the new log-likelihood (depending on convergence criterion) - curr = log_likelihood(matrix, counts) if convergence_criterion=="likelihood" else counts + if convergence_criterion=="likelihood": + curr = log_likelihood(matrix, counts) + else: + curr = counts # Check for convergence - diff = np.abs(curr-prev) if convergence_criterion=="likelihood" else np.abs(curr-prev).sum() + if convergence_criterion=="likelihood": + diff = np.abs(curr-prev) + else: + diff = np.abs(curr-prev).sum() if diff < tolerance: converged = True break diff --git a/irescue/main.py b/irescue/main.py index 8bf3bb1..ceccf3c 100644 --- a/irescue/main.py +++ b/irescue/main.py @@ -186,7 +186,7 @@ def parseArguments(): "--convergence-criterion", type=str, metavar="STR", - choices=["likelihood", "parameters"], + choices=["likelihood", "estimates"], default="likelihood", help=( "Criterion to define convergence. " From ed6974d825ba9780128b5c559a547f056f1996ff Mon Sep 17 00:00:00 2001 From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com> Date: Wed, 15 Apr 2026 14:47:10 +0200 Subject: [PATCH 5/6] update license year --- LICENSE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LICENSE b/LICENSE index 9618698..87748b0 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2022-2025 Benedetto Polimeni +Copyright (c) 2022-2026 Benedetto Polimeni Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal From 60067ea12888ef8aa209d3c151ef598d2ede6e0a Mon Sep 17 00:00:00 2001 From: Benedetto Polimeni <34317613+bepoli@users.noreply.github.com> Date: Wed, 15 Apr 2026 15:29:42 +0200 Subject: [PATCH 6/6] potential bugfix: make copies of counts arrays instead of referencing --- irescue/em.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/irescue/em.py b/irescue/em.py index 9b0187a..d42c859 100644 --- a/irescue/em.py +++ b/irescue/em.py @@ -67,7 +67,7 @@ def run_em(matrix, cycles=100, tolerance=1e-4, convergence_criterion="likelihood if convergence_criterion=="likelihood": prev = log_likelihood(matrix, counts) else: - prev = counts + prev = counts.copy() converged = False curr_cycle = 0 @@ -82,7 +82,7 @@ def run_em(matrix, cycles=100, tolerance=1e-4, convergence_criterion="likelihood if convergence_criterion=="likelihood": curr = log_likelihood(matrix, counts) else: - curr = counts + curr = counts.copy() # Check for convergence if convergence_criterion=="likelihood": @@ -93,6 +93,6 @@ def run_em(matrix, cycles=100, tolerance=1e-4, convergence_criterion="likelihood converged = True break - prev = curr + prev = curr.copy() return counts, (curr_cycle, converged, curr, diff)