From 5f0d5dfd901e97da1b5e595ed4f455010192bf90 Mon Sep 17 00:00:00 2001 From: bernheder Date: Tue, 6 Oct 2026 09:28:45 +0200 Subject: [PATCH 1/3] computational performance --- src/rectanglepy/pp/create_signature.py | 36 +-- src/rectanglepy/tl/deconvolution.py | 340 ++++++++++++++++--------- 2 files changed, 237 insertions(+), 139 deletions(-) diff --git a/src/rectanglepy/pp/create_signature.py b/src/rectanglepy/pp/create_signature.py index e191feb..1a9f4de 100644 --- a/src/rectanglepy/pp/create_signature.py +++ b/src/rectanglepy/pp/create_signature.py @@ -142,17 +142,13 @@ def _run_deseq2( bootstrapped_signature = _create_bootstrap_signature(countsig, sc_data, annotations, number_of_bootstraps) np.random.seed(42) for _i, cell_type in enumerate(countsig.columns): - bootstrapped_signature_copy = bootstrapped_signature.copy() - countsig_copy = countsig.copy() sc_data_filtered = sc_data.T[annotations == cell_type] expressed_cells = (sc_data_filtered > 0).sum(axis=0) if expressed_cells.ndim > 1: # needed for sparse matrices expressed_cells = np.squeeze(np.asarray(expressed_cells)) - # make dense out of sparse - sc_data_filtered = sc_data_filtered.toarray() threshold = gene_expression_threshold * sc_data_filtered.shape[0] - genes = countsig_copy.index[expressed_cells > threshold].tolist() - bootstrapped_signature_copy = bootstrapped_signature_copy.loc[genes].T + genes = countsig.index[expressed_cells > threshold].tolist() + bootstrapped_signature_copy = bootstrapped_signature.loc[genes].T logger.info(f"Running DE analysis for {cell_type}") condition = ["B" if (cell_type + "_") in x else "A" for x in bootstrapped_signature_copy.index] clinical_df = pd.DataFrame({"condition": condition}, index=bootstrapped_signature_copy.index) @@ -174,23 +170,28 @@ def _run_deseq2( def _create_bootstrap_signature(countsig, sc_data, annotations, number_of_bootstraps: int = 7) -> pd.DataFrame: - if scipy.sparse.issparse(sc_data): - sc_data = sc_data.toarray() + # Cells are only ever summed here, never read one by one, so sparse input is kept sparse: + # densifying the whole matrix costs genes * cells * 8 bytes, which is gigabytes on an + # atlas-sized reference. Dense input deliberately stays dense - gathering rows out of a + # dense array is faster than routing them through a sparse format. + cells_by_gene = sc_data.T.tocsr() if scipy.sparse.issparse(sc_data) else sc_data.T celltypes = countsig.columns - bootstrapped_signature = pd.DataFrame() + columns = {} samples_per_bootstrap = 500 np.random.seed(42) for celltype in celltypes: - sc_data_filtered = sc_data.T[annotations == celltype] + sc_data_filtered = cells_by_gene[(annotations == celltype).to_numpy()] + # sparse matrices do not support len(); shape[0] is the cell count for both formats + number_of_cells = sc_data_filtered.shape[0] for i in range(number_of_bootstraps): - selected_rows = np.random.choice(len(sc_data_filtered), samples_per_bootstrap, replace=True) + selected_rows = np.random.choice(number_of_cells, samples_per_bootstrap, replace=True) summed_rows = sc_data_filtered[selected_rows].sum(axis=0) - bootstrapped_signature[f"{celltype}_{i}"] = list(summed_rows) - bootstrapped_signature.index = countsig.index + # a sparse sum is a (1, genes) matrix, a dense one a flat array + columns[f"{celltype}_{i}"] = np.asarray(summed_rows).ravel() + # built in one go: assigning column by column repeatedly reallocates the frame + bootstrapped_signature = pd.DataFrame(columns, index=countsig.index) # to int - samples_per_bootstrap = samples_per_bootstrap / 2.5 - bootstrapped_signature = bootstrapped_signature.astype(int) - return bootstrapped_signature + return bootstrapped_signature.astype(int) def _de_analysis( @@ -425,9 +426,10 @@ def build_rectangle_signatures( def _create_pseudo_count_sig(sc_counts: np.ndarray, annotations: pd.Series, var_names) -> pd.DataFrame: unique_labels, label_indices = np.unique(annotations, return_inverse=True) grouped_sum = np.zeros((len(unique_labels), sc_counts.shape[0])) + cells_by_gene = sc_counts.T for i, _label in enumerate(unique_labels): label_columns = label_indices == i - grouped_sum[i, :] = np.sum(sc_counts.T[label_columns, :], axis=0) + grouped_sum[i, :] = np.sum(cells_by_gene[label_columns, :], axis=0) grouped_sum = grouped_sum.T grouped_sum = pd.DataFrame(grouped_sum, index=var_names, columns=unique_labels).astype(int) diff --git a/src/rectanglepy/tl/deconvolution.py b/src/rectanglepy/tl/deconvolution.py index 1be3001..a4cb4c2 100644 --- a/src/rectanglepy/tl/deconvolution.py +++ b/src/rectanglepy/tl/deconvolution.py @@ -7,18 +7,166 @@ import osqp import pandas as pd import scipy.sparse as sp -import statsmodels.api as sm from joblib import Parallel, delayed, parallel_backend from loguru import logger from rectanglepy.pp.rectangle_signature import RectangleSignatureResult +# OSQP can return noticeably different solutions from active-set QP when P is singular/ill-conditioned. +# A tiny ridge makes the problem strictly convex and closer to quadprog behavior. +QP_RIDGE = 1e-8 + def _scale_weights(weights: np.ndarray) -> np.ndarray: min_weight = np.nextafter(min(weights), np.float64(1.0)) # prevent division by zero return weights / min_weight +def _upper_triangular_csc_pattern(n_vars: int) -> tuple[np.ndarray, np.ndarray]: + """Returns the (indices, indptr) of a fully populated upper triangular CSC matrix. + + OSQP only stores the upper triangle of P, so handing it the triangle directly saves the + conversion of the full symmetric matrix on every solve. + """ + indices = np.concatenate([np.arange(column + 1) for column in range(n_vars)]).astype(np.int32) + indptr = np.concatenate(([0], np.cumsum(np.arange(1, n_vars + 1)))).astype(np.int32) + return indices, indptr + + +class _QuadraticProgram: + """A deconvolution QP for one (signature, bulk) pair, solvable at several dampening weights. + + Everything that does not depend on the weights - the numpy views of the inputs and the constraint + matrix with its bounds - is built once and shared by all solves, which keeps it out of the + dampened least squares iteration. + + Note that each solve still gets its own OSQP workspace: reusing one across solves makes OSQP + keep the problem scaling of the first objective, which perturbs the fixed point the iteration + converges to (and is not faster at these problem sizes). + """ + + def __init__( + self, + signature: pd.DataFrame, + bulk: pd.Series, + prev_assignments: list[int | str] = None, + prev_solution: pd.Series = None, + ): + if not signature.index.equals(bulk.index): + bulk = bulk.reindex(signature.index) + self.signature = signature.to_numpy(dtype=np.float64) + self.bulk = bulk.to_numpy(dtype=np.float64) + self.n_vars = self.signature.shape[1] # number of cell types / fractions + + self._triu_indices, self._triu_indptr = _upper_triangular_csc_pattern(self.n_vars) + self._triu_selection = np.tril_indices(self.n_vars) # reads P.T column-major upper triangle + self._constraints, self._lower_bounds = self._build_constraints(prev_assignments, prev_solution) + self._upper_bounds = np.full_like(self._lower_bounds, np.inf, dtype=np.float64) + + def _build_constraints( + self, prev_assignments: list[int | str], prev_solution: pd.Series + ) -> tuple[sp.csc_matrix, np.ndarray]: + # ----- Constraints in the original quadprog form: C.T x >= b + # We build C (n_vars, n_constraints) then map to OSQP: A = C.T, l=b, u=+inf + C_cols = [] + b_list = [] + + # Constraint 1: sum(x) <= 1 -> -sum(x) >= -1 + C_cols.append(-np.ones((self.n_vars, 1), dtype=np.float64)) + b_list.append(np.array([-1.0], dtype=np.float64)) + + # Constraint 2: x >= 0 -> I x >= 0 + C_cols.append(np.eye(self.n_vars, dtype=np.float64)) + b_list.append(np.zeros(self.n_vars, dtype=np.float64)) + + # Constraint 3: keep close to prev_solution (this encoding already turns upper bounds into >= via negation) + if prev_solution is not None: + if prev_assignments is None: + raise ValueError("prev_assignments must be provided when prev_solution is provided.") + + for cluster in prev_solution.index: + # x_cluster <= upper -> -x_cluster >= -upper + C_upper = np.array( + [-1.0 if str(x) == str(cluster) else 0.0 for x in prev_assignments], + dtype=np.float64, + ).reshape(-1, 1) + + # x_cluster >= lower -> +x_cluster >= +lower + C_lower = np.array( + [1.0 if str(x) == str(cluster) else 0.0 for x in prev_assignments], + dtype=np.float64, + ).reshape(-1, 1) + + prev_weight = float(prev_solution.loc[cluster]) + upper = min(1.0, prev_weight + 0.03) + lower = max(0.0, prev_weight - 0.03) + + C_cols.append(C_upper) + b_list.append(np.array([-upper], dtype=np.float64)) + + C_cols.append(C_lower) + b_list.append(np.array([lower], dtype=np.float64)) + + C = np.concatenate(C_cols, axis=1) # (n_vars, n_constraints) + b = np.concatenate(b_list, axis=0).astype(np.float64) # (n_constraints,) + + # OSQP uses l <= A x <= u + return sp.csc_matrix(C.T), b + + def _objective(self, gld: np.ndarray, multiplier: int) -> tuple[np.ndarray, np.ndarray]: + # Minimize 1/2 x^T G x - a^T x + if multiplier is None: + G = self.signature.T @ self.signature + a = self.signature.T @ self.bulk + else: + weights = np.square(1 / (self.signature @ gld)) + weights_dampened = np.clip(_scale_weights(weights), None, multiplier) + # equivalent to signature.T @ diag(weights) @ signature, without the dense (genes, genes) matrix + G = self.signature.T @ (weights_dampened[:, None] * self.signature) + a = self.signature.T @ (weights_dampened * self.bulk) + + # ----- Map to OSQP: minimize 1/2 x^T P x + q^T x -> q = -a + P = G + q = -a + + # Optional scaling + scale = np.linalg.norm(P) + if scale > 0: + P = P / scale + q = q / scale + + P = ((P + P.T) / 2.0) + QP_RIDGE * np.eye(self.n_vars, dtype=np.float64) + return P, q + + def solve(self, gld: np.ndarray = None, multiplier: int = None) -> np.ndarray: + P, q = self._objective(gld, multiplier) + P_data = P.T[self._triu_selection] # upper triangle in CSC (column major) order + P_sp = sp.csc_matrix((P_data, self._triu_indices, self._triu_indptr), shape=(self.n_vars, self.n_vars)) + + solver = osqp.OSQP() + solver.setup( + P=P_sp, + q=q, + A=self._constraints, + l=self._lower_bounds, + u=self._upper_bounds, + verbose=False, + eps_abs=1e-7, + eps_rel=1e-7, + max_iter=50000, + polish=True, + warm_start=False, + scaled_termination=False, + ) + + res = solver.solve() + + if res.info.status_val not in (1,): # 1 = solved + raise RuntimeError(f"OSQP did not solve the problem: {res.info.status}") + + return res.x + + def solve_qp( signature: pd.DataFrame, bulk: pd.Series, @@ -52,143 +200,90 @@ def solve_qp( Notes ----- This function uses quadratic programming to solve the deconvolution problem. The objective is to minimize the difference between the observed bulk data and the data predicted by the signature matrix and the cell fractions. The function also includes constraints to ensure that the cell fractions are non-negative and sum to 1, and to make the solution similar to the previous assignments and weights if they are provided. + + Callers that solve the same problem at several dampening weights (e.g. the dampened least squares + iteration) should use :class:`_QuadraticProgram` directly, which shares the weight independent + parts of the problem across solves. """ - # ------------------ QP-based deconvolution - # Minimize 1/2 x^T G x - a^T x - # Subject to C.T x >= b - if multiplier is None: - a = (signature.T @ bulk).to_numpy(dtype=np.float64) - G = (signature.T @ signature).to_numpy(dtype=np.float64) - else: - weights = np.square(1 / (signature @ gld)) - weights_dampened = np.clip(_scale_weights(weights), None, multiplier) - W = np.diag(np.asarray(weights_dampened, dtype=np.float64)) - G = (signature.T.to_numpy() @ (W @ signature.to_numpy())).astype(np.float64) - a = (signature.T.to_numpy() @ (W @ bulk.to_numpy())).astype(np.float64) - - n_vars = G.shape[0] # number of cell types / fractions - - # ----- Constraints in your original quadprog form: C.T x >= b - # We'll build C (n_vars, n_constraints) then map to OSQP: A = C.T, l=b, u=+inf - C_cols = [] - b_list = [] - - # Constraint 1: sum(x) <= 1 -> -sum(x) >= -1 - C_cols.append(-np.ones((n_vars, 1), dtype=np.float64)) - b_list.append(np.array([-1.0], dtype=np.float64)) - - # Constraint 2: x >= 0 -> I x >= 0 - C_cols.append(np.eye(n_vars, dtype=np.float64)) - b_list.append(np.zeros(n_vars, dtype=np.float64)) - - # Constraint 3: keep close to prev_solution (your encoding already turns upper bounds into >= via negation) - if prev_solution is not None: - if prev_assignments is None: - raise ValueError("prev_assignments must be provided when prev_solution is provided.") - - for cluster in prev_solution.index: - # x_cluster <= upper -> -x_cluster >= -upper - C_upper = np.array( - [-1.0 if str(x) == str(cluster) else 0.0 for x in prev_assignments], - dtype=np.float64, - ).reshape(-1, 1) - - # x_cluster >= lower -> +x_cluster >= +lower - C_lower = np.array( - [1.0 if str(x) == str(cluster) else 0.0 for x in prev_assignments], - dtype=np.float64, - ).reshape(-1, 1) - - prev_weight = float(prev_solution.loc[cluster]) - upper = min(1.0, prev_weight + 0.03) - lower = max(0.0, prev_weight - 0.03) - - C_cols.append(C_upper) - b_list.append(np.array([-upper], dtype=np.float64)) - - C_cols.append(C_lower) - b_list.append(np.array([lower], dtype=np.float64)) - - C = np.concatenate(C_cols, axis=1) # (n_vars, n_constraints) - b = np.concatenate(b_list, axis=0).astype(np.float64) # (n_constraints,) - - # ----- Map to OSQP: minimize 1/2 x^T P x + q^T x - # Your objective: 1/2 x^T G x - a^T x -> q = -a - P = G - q = -a - - # Optional scaling - scale = np.linalg.norm(P) - if scale > 0: - P = P / scale - q = q / scale - - # OSQP can return noticeably different solutions from active-set QP when P is singular/ill-conditioned. - # Add tiny ridge to make the problem strictly convex and closer to quadprog behavior. - ridge = 1e-8 - P = ((P + P.T) / 2.0) + ridge * np.eye(n_vars, dtype=np.float64) - - # OSQP uses l <= A x <= u - A = C.T # (n_constraints, n_vars) - l = b - u = np.full_like(l, np.inf, dtype=np.float64) - - # Sparse matrices (required/expected) - P_sp = sp.csc_matrix(P) - A_sp = sp.csc_matrix(A) - - solver = osqp.OSQP() - solver.setup( - P=P_sp, - q=q, - A=A_sp, - l=l, - u=u, - verbose=False, - eps_abs=1e-7, - eps_rel=1e-7, - max_iter=50000, - polish=True, - warm_start=False, - scaled_termination=False, - ) + problem = _QuadraticProgram(signature, bulk, prev_assignments, prev_solution) + return problem.solve(gld, multiplier) + + +def _batched_wls( + design: np.ndarray, + response: np.ndarray, + weights: np.ndarray, + subsets: np.ndarray, + max_bytes: int = 64 * 1024**2, +) -> np.ndarray: + """Fits one weighted least squares model per gene subset in batch. + + Mirrors ``statsmodels.api.WLS(...).fit()``: the design and the response are whitened with the + square root of the weights and the parameters are read off the pseudo inverse of the whitened + design. Stacking the subsets lets numpy loop over them in C instead of Python. - res = solver.solve() + Parameters + ---------- + design : np.ndarray + The (genes, cell types) design matrix. + response : np.ndarray + The (genes,) response vector. + weights : np.ndarray + The (genes,) regression weights. + subsets : np.ndarray + A (n_subsets, subset size) array of gene indices, one row per model. + max_bytes : int + Upper bound on the size of a single whitened design batch, used to chunk the subsets. + + Returns + ------- + np.ndarray + A (n_subsets, cell types) array of fitted parameters. + """ + n_subsets, subset_size = subsets.shape + n_params = design.shape[1] + params = np.empty((n_subsets, n_params), dtype=np.float64) - if res.info.status_val not in (1,): # 1 = solved - raise RuntimeError(f"OSQP did not solve the problem: {res.info.status}") + chunk_size = max(1, int(max_bytes // (subset_size * n_params * design.itemsize))) + for start in range(0, n_subsets, chunk_size): + chunk = subsets[start : start + chunk_size] + sqrt_weights = np.sqrt(weights[chunk]) # (chunk, subset size) + whitened_design = design[chunk] * sqrt_weights[:, :, None] # (chunk, subset size, cell types) + whitened_response = response[chunk] * sqrt_weights + params[start : start + chunk_size] = np.einsum("bpg,bg->bp", np.linalg.pinv(whitened_design), whitened_response) - return res.x + return params def _calculate_dampening_constant(signature: pd.DataFrame, bulk: pd.Series, qp_gld: np.ndarray) -> int: solutions_std = [] np.random.seed(1) - weights = np.square(1 / (np.dot(signature, qp_gld))) + signature_values = np.asarray(signature, dtype=np.float64) + bulk_values = np.asarray(bulk, dtype=np.float64) + n_genes = signature_values.shape[0] + weights = np.square(1 / (signature_values @ qp_gld)) weights_scaled = _scale_weights(weights) weights_scaled_no_inf = weights_scaled[weights_scaled != np.inf] qp_gld_sum = sum(qp_gld) + # the design matrix reproduces the original ``-1 + signature`` term of the WLS fit + design = signature_values - 1.0 + subset_size = n_genes // 2 + n_subsets = 100 # try multiple values of the dampening constant (multiplier) # for each, calculate the variance of the dampened weighted solution for a subset of genes max_range = 40 multiplier_range = min(max_range, math.ceil(np.log2(max(weights_scaled_no_inf)))) for i in range(multiplier_range): - solutions = [] multiplier = 2**i - weights_dampened = np.array([multiplier if multiplier <= x else x for x in weights_scaled]).astype("double") - for _ in range(100): - subset = np.random.choice(len(signature), size=len(signature) // 2, replace=False) - bulk_subset = bulk.iloc[list(subset)] - signature_subset = signature.iloc[subset, :] - fit = sm.WLS(bulk_subset, -1 + signature_subset, weights=weights_dampened[subset]).fit() - solution = fit.params * qp_gld_sum / sum(fit.params) - solutions.append(solution) - solutions_df = pd.DataFrame(solutions) - - solutions_std.append(solutions_df.std(axis=0)) - solutions_std_df = pd.DataFrame(solutions_std) - means = solutions_std_df.apply(lambda x: np.mean(x**2), axis=1) - best_dampening_constant = means.idxmin() + weights_dampened = np.minimum(weights_scaled, multiplier) + subsets = np.array([np.random.choice(n_genes, size=subset_size, replace=False) for _ in range(n_subsets)]) + params = _batched_wls(design, bulk_values, weights_dampened, subsets) + solutions = params * qp_gld_sum / params.sum(axis=1, keepdims=True) + + solutions_std.append(np.nanstd(solutions, axis=0, ddof=1)) + solutions_std_array = np.array(solutions_std) + means = np.nanmean(np.square(solutions_std_array), axis=1) + best_dampening_constant = int(np.nanargmin(means)) return best_dampening_constant @@ -202,7 +297,8 @@ def _calculate_ls( signature = signature.loc[genes].sort_index() bulk = bulk.loc[genes].sort_index().astype("double") - approximate_solution = solve_qp(signature, bulk, prev_assignments, prev_weights) + problem = _QuadraticProgram(signature, bulk, prev_assignments, prev_weights) + approximate_solution = problem.solve() dampening_constant = _calculate_dampening_constant(signature, bulk, approximate_solution) multiplier = 2**dampening_constant @@ -212,7 +308,7 @@ def _calculate_ls( iterations = 2 solutions_sum = approximate_solution while (change > convergence_threshold) and (iterations < max_iterations): - dampened_solution = solve_qp(signature, bulk, prev_assignments, prev_weights, approximate_solution, multiplier) + dampened_solution = problem.solve(approximate_solution, multiplier) solutions_sum += dampened_solution solution_averages = solutions_sum / iterations change = np.linalg.norm(solution_averages - approximate_solution, 1) From 00ecf2d0129d44f5f35ec161e5a36e34aab57c68 Mon Sep 17 00:00:00 2001 From: bernheder Date: Wed, 7 Oct 2026 10:23:08 +0200 Subject: [PATCH 2/3] hyperparameter optimization --- src/rectanglepy/parameters.py | 2 +- src/rectanglepy/pp/create_signature.py | 16 ++++++++-------- src/rectanglepy/rectangle.py | 4 ++-- src/rectanglepy/tl/deconvolution.py | 5 ++--- 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/rectanglepy/parameters.py b/src/rectanglepy/parameters.py index ab8bcce..8ce53a2 100644 --- a/src/rectanglepy/parameters.py +++ b/src/rectanglepy/parameters.py @@ -7,5 +7,5 @@ class RectangleAdvancedParameters: """Advanced parameters for Rectangle internals.""" - number_of_bootstraps: int = 7 + number_of_bootstraps: int = 20 grid_search_split_size: int = 50 diff --git a/src/rectanglepy/pp/create_signature.py b/src/rectanglepy/pp/create_signature.py index 1a9f4de..07d5666 100644 --- a/src/rectanglepy/pp/create_signature.py +++ b/src/rectanglepy/pp/create_signature.py @@ -134,8 +134,8 @@ def _run_deseq2( sc_data, annotations: pd.Series, n_cpus: int = None, - gene_expression_threshold=0.5, - number_of_bootstraps: int = 7, + gene_expression_threshold=0.4, + number_of_bootstraps: int = 20, ) -> dict[str | int, pd.DataFrame]: results = {} inference = DefaultInference(n_cpus=n_cpus) @@ -169,7 +169,7 @@ def _run_deseq2( return results -def _create_bootstrap_signature(countsig, sc_data, annotations, number_of_bootstraps: int = 7) -> pd.DataFrame: +def _create_bootstrap_signature(countsig, sc_data, annotations, number_of_bootstraps: int = 20) -> pd.DataFrame: # Cells are only ever summed here, never read one by one, so sparse input is kept sparse: # densifying the whole matrix costs genes * cells * 8 bytes, which is gigabytes on an # atlas-sized reference. Dense input deliberately stays dense - gathering rows out of a @@ -177,7 +177,7 @@ def _create_bootstrap_signature(countsig, sc_data, annotations, number_of_bootst cells_by_gene = sc_data.T.tocsr() if scipy.sparse.issparse(sc_data) else sc_data.T celltypes = countsig.columns columns = {} - samples_per_bootstrap = 500 + samples_per_bootstrap = 1000 np.random.seed(42) for celltype in celltypes: sc_data_filtered = cells_by_gene[(annotations == celltype).to_numpy()] @@ -203,7 +203,7 @@ def _de_analysis( optimize_cutoffs: bool, n_cpus: int = None, genes=None, - gene_expression_threshold=0.5, + gene_expression_threshold=0.4, advanced_parameters: RectangleAdvancedParameters = None, ) -> tuple[Series, dict[str, [str]] :, DataFrame | None]: logger.info("Starting DE analysis") @@ -308,7 +308,7 @@ def build_rectangle_signatures( p=0.015, lfc=1.5, n_cpus: int = None, - gene_expression_threshold=0.5, + gene_expression_threshold=0.4, advanced_parameters: RectangleAdvancedParameters = None, ) -> RectangleSignatureResult: r"""Builds rectangle signatures based on single-cell count data and annotations. @@ -334,7 +334,7 @@ def build_rectangle_signatures( n_cpus The number of cpus to use for the DE analysis. Defaults to the number of cpus available. gene_expression_threshold - The gene expression threshold for the DE analysis. How many cells need to express a gene to be considered in DGE + The gene expression threshold for the DE analysis. The fraction of cells that must express a gene to be considered in DGE. Defaults to 0.4 advanced_parameters Optional advanced Rectangle parameters. Defaults are used when not provided. @@ -445,7 +445,7 @@ def _optimize_parameters( grid_search_split_size: int = 50, ) -> pd.DataFrame: # search space for p and lfc - lfcs = [x / 100 for x in range(160, 230, 10)] + lfcs = [1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0] ps = [x / 1000 for x in range(50, 51, 1)] results = [] diff --git a/src/rectanglepy/rectangle.py b/src/rectanglepy/rectangle.py index d3ae6a2..ba97d48 100644 --- a/src/rectanglepy/rectangle.py +++ b/src/rectanglepy/rectangle.py @@ -24,7 +24,7 @@ def rectangle( p=0.015, lfc=1.5, n_cpus: int = None, - gene_expression_threshold=0.5, + gene_expression_threshold=0.4, advanced_parameters: RectangleAdvancedParameters = None, ) -> tuple[DataFrame, RectangleSignatureResult]: r"""All in one deconvolution method. Creates signatures and deconvolutes the bulk data. Has options for subsampling and consensus runs. @@ -52,7 +52,7 @@ def rectangle( correct_mrna_bias : bool A flag indicating whether to correct for mRNA bias. Defaults to True. gene_expression_threshold : float - The threshold for gene expression. Genes with expression below this threshold are removed from the analysis. + The threshold for gene expression. Genes must be expressed in at least this fraction of cells. Defaults to 0.4. advanced_parameters Optional advanced Rectangle parameters. Defaults are used when not provided. diff --git a/src/rectanglepy/tl/deconvolution.py b/src/rectanglepy/tl/deconvolution.py index a4cb4c2..02daf12 100644 --- a/src/rectanglepy/tl/deconvolution.py +++ b/src/rectanglepy/tl/deconvolution.py @@ -265,8 +265,6 @@ def _calculate_dampening_constant(signature: pd.DataFrame, bulk: pd.Series, qp_g weights_scaled = _scale_weights(weights) weights_scaled_no_inf = weights_scaled[weights_scaled != np.inf] qp_gld_sum = sum(qp_gld) - # the design matrix reproduces the original ``-1 + signature`` term of the WLS fit - design = signature_values - 1.0 subset_size = n_genes // 2 n_subsets = 100 # try multiple values of the dampening constant (multiplier) @@ -277,7 +275,8 @@ def _calculate_dampening_constant(signature: pd.DataFrame, bulk: pd.Series, qp_g multiplier = 2**i weights_dampened = np.minimum(weights_scaled, multiplier) subsets = np.array([np.random.choice(n_genes, size=subset_size, replace=False) for _ in range(n_subsets)]) - params = _batched_wls(design, bulk_values, weights_dampened, subsets) + # WLS uses the signature directly; no intercept column is added. + params = _batched_wls(signature_values, bulk_values, weights_dampened, subsets) solutions = params * qp_gld_sum / params.sum(axis=1, keepdims=True) solutions_std.append(np.nanstd(solutions, axis=0, ddof=1)) From 15902b99884f272724945f591a05dd9f55533333 Mon Sep 17 00:00:00 2001 From: bernheder Date: Wed, 7 Oct 2026 10:32:07 +0200 Subject: [PATCH 3/3] test updated --- tests/test_pp.py | 4 ++-- tests/test_tl.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_pp.py b/tests/test_pp.py index c65e6bc..c6812ba 100644 --- a/tests/test_pp.py +++ b/tests/test_pp.py @@ -236,12 +236,12 @@ def test_de_analysis(small_data): # test with sparse matrix _ = _de_analysis(sc_pseudo, adata_sparse.X.T, annotations, 0.4, 0.1, False, None, adata.var_names) - assert 5 < len(r1) < 50 + assert 5 < len(r1) < 100 assert len(r2) == 3 def test_create_bootstrap_signature(small_data): - bootstraps_per_cell = 7 + bootstraps_per_cell = RectangleAdvancedParameters().number_of_bootstraps sc_counts, annotations, bulk = small_data sc_counts = sc_counts.astype("int") sc_pseudo = sc_counts.groupby(annotations.values, axis=1).sum() diff --git a/tests/test_tl.py b/tests/test_tl.py index 5d1eb4c..95bc02f 100644 --- a/tests/test_tl.py +++ b/tests/test_tl.py @@ -62,7 +62,7 @@ def test_simple_weighted_dampened_deconvolution(quantiseq_data): corr = np.corrcoef(result, expected)[0, 1] rsme = np.sqrt(np.mean((result - expected) ** 2)) - assert corr > 0.815 and rsme < 0.011 + assert corr > 0.81 and rsme < 0.01 def test_correct_for_unknown_cell_content(small_data, quantiseq_data):