From d9902f31fd725cc237efc4461b29fa024ee66c5e Mon Sep 17 00:00:00 2001 From: Marco Cominelli <161143279+marco-cominelli01@users.noreply.github.com> Date: Wed, 17 Dec 2025 21:43:34 +0100 Subject: [PATCH 1/3] Add noisy-feature handling and extra convergence criterion --- irescue/main.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/irescue/main.py b/irescue/main.py index dc7d5e5..5b33fac 100644 --- a/irescue/main.py +++ b/irescue/main.py @@ -178,6 +178,25 @@ def parseArguments(): "(Default: %(default)s)." ), ) + parser.add_argument( + "--keep-all-features", + action="store_true", + help=( + "Perform EM considering 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", @@ -376,6 +395,8 @@ def main(): feature_index, dirs["tmp"], args.no_umi, + args.keep_all_features, + args.convergence_criterion, args.dump_ec, args.max_iters, args.tolerance, From d44e545518705858fec8d668ce402e08a153d836 Mon Sep 17 00:00:00 2001 From: Marco Cominelli <161143279+marco-cominelli01@users.noreply.github.com> Date: Thu, 18 Dec 2025 15:22:42 +0100 Subject: [PATCH 2/3] Enhance compute_cell_counts to handle noisy loci and make available another convergence criterion --- irescue/count.py | 42 ++++++++++++++++++++++++++++++------------ 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/irescue/count.py b/irescue/count.py index 2f544f2..480a240 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, keep_all_features, convergence_criterion ): """ Calculate TE counts of a single cell, given a list of equivalence classes. @@ -259,18 +259,32 @@ 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 keep_all_features: + tokeep = np.flatnonzero(em_array.sum(axis=0)) + else: + # keep features supported by at least 2 multimapping reads + tokeep1 = np.where((em_array.sum(axis=0) >= 2).A1)[0] + # rescue features supported by only 1 multimapping read but also by (at least) 1 uniquely mapping read + 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 step if keep_all_features=True) + if not keep_all_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 +306,8 @@ def run_count( features_index, tmpdir, no_umi, + keep_all_features, + convergence_criterion, dumpEC, max_iters, tolerance, @@ -323,11 +339,13 @@ def run_count( tolerance=tolerance, dumpEC=dumpEC, no_umi=no_umi, + keep_all_features=keep_all_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, ) From 4af6723b0d83db25f6e217b2122582f170c44c02 Mon Sep 17 00:00:00 2001 From: Marco Cominelli <161143279+marco-cominelli01@users.noreply.github.com> Date: Thu, 18 Dec 2025 15:28:41 +0100 Subject: [PATCH 3/3] Enhance run_em with convergence criterion parameter --- irescue/em.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/irescue/em.py b/irescue/em.py index 32c94e7..f451522 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 parameters) + 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 (or new parameters) + curr = log_likelihood(matrix, counts) if convergence_criterion=="likelihoods" 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)