Expanded ebayes ttest by vbrennsteiner · Pull Request #240 · MannLabs/alphapepttools

vbrennsteiner · 2026-06-26T20:52:57Z

Purpose of this PR

The existing implementation of diff_exp_ebayes() works well for differential expression data but has some critical shortcomings, the first of which is directly caused by the underlying lmFit() implementation of inmoose:

No handling of NaN-containing data, i.e. requiring 100 % feature completeness to avoid crashing lmFit()
Only computing one contrast (e.g. Compound_X vs. Control) at a time
Not supporting any linear covariates that could be added to the linear fit to explain batch variance

ebayes_expanded.py addresses all of these points in its new version of diff_exp_ebayes, which consists of six new custom functions to sidestep all inmoose dependencies except limma.ebayes itself. Its critical component is nan_lmfit(), performs the same kind of linear fitting as lmFit() but tolerates missing values:

build_design_matrix to construct a 1/0 matrix for fitting the linear coefficients
generate_feature_mask to construct a selection mask for features with sufficient control samples (some tests experimental setups require an up-front check to only assess features with sufficient control replicates. This function separates that logic from the actual linear fitting below
nan_lmfit to perform the actual feature-wise fitting while tolerating missing features. Can optionally take a feature_mask from generate_feature_mask
make_contrasts constructs a contrasts matrix to construct fold changes from the coefficient's delta (i.e. Beta_treatment - Beta_control) results in the respective fold-change
run_contrasts Actually compute the fold changes, i.e. the t-test but keep unscaled variance and standard deviation in a separate vector; Since the coefficients (i.e. also the coefficient's delta's) standard error consists of the standard deviation and the variance, and empirical Bayes moderation adjusts the standard deviation, these two must remain separate here
ebayes_moderation takes the coefficients and the unscaled variance and standard deviation arrays and performs empirical Bayes moderation. This results in modified t-values and p-values, and is also the only point where we still need inmoose

These scripts are coordinated in the new version of diff_exp_ebayes, which can now handle nans, multiple contrasts and linear covariates, i.e. to model batch effects.

Critical tests

The key test is numerical equivalence between the new and old diff_exp_ebayes versions in test_diff_exp.py on a vanilla dataset without missing values or covariates and with only one contrast. Furthermore, all new functions are covered by tests, with special attention on nan_lmfit to ensure that there are no index misalignments of features or coefficients when arrays are brought back to full shape after fitting with filtered design matrices that avoid missing samples for the respective feature

Outlook

The new function can be used largely as-is in place of the existing diff_exp_ebayes with one breaking change: results are now returned in a {contrast_name: dataframe} structured dict, rather than a simple dataframe in the previous single-comparison case. This is a direct result of supporting multiple contrasts now.

Once a port of the limma.ebayes function from inmoose exists, inmoose can be dropped as a dependency entirely.

…test

…sting implementation#

… special focus on preventing any misalignments during the nan_lfit step, wehre the matrices with skipped coefficients are scattered back to full size

…d on control samples, and nan_lmfit that purely handles linear fitting on the adata with a design matrix and an optional feature mask

…s gated with MNAR, where even if there are too few samples to report a fold change, the samples that are there are still used by the Bayesian model. however, the high condition is gated by high_max_missing, which excludes features altogether if their high-signal condition has too many missing values. The logic here is that if a feature drops out randomly when the signal should be high, it is unreliable and we don't want to use it for modeling

…reatment

codecov-commenter · 2026-06-29T00:49:05Z

Codecov Report

❌ Patch coverage is 86.27451% with 28 lines in your changes missing coverage. Please review.
✅ Project coverage is 73.94%. Comparing base (ca7b133) to head (a09334e).
⚠️ Report is 118 commits behind head on main.

Files with missing lines	Patch %	Lines
src/alphapepttools/tl/diff_exp/ebayes_expanded.py	86.27%	28 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #240      +/-   ##
==========================================
+ Coverage   73.19%   73.94%   +0.74%     
==========================================
  Files          35       36       +1     
  Lines        2067     2268     +201     
==========================================
+ Hits         1513     1677     +164     
- Misses        554      591      +37

Files with missing lines	Coverage Δ
src/alphapepttools/tl/utils.py	`85.71% <ø> (-0.57%)`	⬇️
src/alphapepttools/tl/diff_exp/ebayes_expanded.py	`86.27% <86.27%> (ø)`

... and 5 files with indirect coverage changes

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

…drop the asymmetrical gates so that differential expression doesn't require knowledge of which side would/could be mnar and which not --> this type of decision is now made upstream with pp.filter_data_completeness for MAR and internally with symmetrical thresholds for MNAR for either comparison side

…me functionality to small helper functions to make filter_data_completeness less complex

…roupby operation

…and cover with unit tests

mschwoer · 2026-07-13T14:25:51Z

        )


+def _validate_max_missing(


this is #242, right?

mschwoer · 2026-07-13T14:26:19Z

+)
+
+
+def build_design_matrix(


could this be private? (check also the rest of the public API of this module)

mschwoer · 2026-07-13T14:27:44Z

+    return dm, {"condition_col_idxs": condition_col_idxs, "covariate_col_idxs": covariate_col_idxs}
+
+
+def summarize_design_matrix(


print_design_matrix_summary() ?

mschwoer · 2026-07-13T14:34:14Z

+    if between_column not in adata.obs.columns:
+        raise ValueError(f"Column '{between_column}' not found in adata.obs.")
+    between_levels = adata.obs[between_column].unique()
+
+    # Validate the B condition (the single reference, comparison[1])
+    b_condition = comparison[1]
+    if b_condition not in between_levels:
+        raise ValueError(f"Condition '{b_condition}' not found in column '{between_column}'.")
+
+    # Validate the A conditions (comparison[0])
+    a_conditions = comparison[0]
+    if a_conditions == "_ALL_":
+        a_conditions = [level for level in between_levels if level != b_condition]
+    elif isinstance(a_conditions, str):
+        a_conditions = [a_conditions]
+
+    for a_condition in a_conditions:
+        if a_condition not in between_levels:
+            raise ValueError(f"Condition '{a_condition}' not found in column '{between_column}'.")


this method is quite long, move this to a dedicated method _validate_input() or smth?

mschwoer · 2026-07-13T14:35:44Z

+    contrast_names = contrasts_from_matrix(contrast_matrix, b_condition)
+    if len(contrast_names) != contrast_results["log2fc"].shape[0]:
+        raise ValueError("Number of contrast names does not match number of contrasts in results.")
+
+    results = {}
+    for contrast_idx, contrast_name in enumerate(contrast_names):
+        # By convention the contrast is named "A_VS_B" (comparison[0] vs comparison[1]).
+        a_level, b_level = contrast_name.split("_VS_")
+
+        p_values = ebayes_results["p"][contrast_idx].copy()
+        log2fc = contrast_results["log2fc"][contrast_idx].copy()
+
+        # Replicate gate: suppress the fold change unless both conditions have enough observed values
+        keep = np.ones(adata.n_vars, dtype=bool)
+        if a_min_required is not None:
+            keep &= _sufficient_values_mask(adata, between_column, a_level, a_min_required)
+        if b_min_required is not None:
+            keep &= _sufficient_values_mask(adata, between_column, b_level, b_min_required)
+        p_values[~keep] = np.nan
+        log2fc[~keep] = np.nan
+
+        # preprocess p-values and FDR for the current contrast
+        fdr_pvalues = nan_safe_bh_correction(p_values)
+        neg_log10_fdr = np.array([negative_log10_pvalue(fdr) for fdr in fdr_pvalues])
+        neg_log10_pvalues = np.array([negative_log10_pvalue(p) for p in p_values])
+
+        # Sample counts per level, mirroring ebayes.py. The name is "A_VS_B".
+        max_level_1_samples, max_level_2_samples = determine_max_replicates(adata, between_column, a_level, b_level)
+
+        results[contrast_name] = pd.DataFrame(
+            {
+                "condition_pair": contrast_name,
+                "protein": adata.var_names,
+                "log2fc": log2fc,
+                "p_value": p_values,
+                "-log10(p_value)": neg_log10_pvalues,
+                "fdr": fdr_pvalues,
+                "-log10(fdr)": neg_log10_fdr,
+                "method": "limma_ebayes_inmoose_expanded",
+                "max_level_1_samples": max_level_1_samples,
+                "max_level_2_samples": max_level_2_samples,
+            },
+            index=adata.var_names,
+        )


this has a different level of abstraction than the rest of this method's code, move it to a dedicated method

mschwoer · 2026-07-13T14:36:27Z

+        log2fc = contrast_results["log2fc"][contrast_idx].copy()
+
+        # Replicate gate: suppress the fold change unless both conditions have enough observed values
+        keep = np.ones(adata.n_vars, dtype=bool)


(nit) keep_mask

mschwoer · 2026-07-13T14:37:14Z

+    results = {}
+    for contrast_idx, contrast_name in enumerate(contrast_names):
+        # By convention the contrast is named "A_VS_B" (comparison[0] vs comparison[1]).
+        a_level, b_level = contrast_name.split("_VS_")


introduce a string constant for _VS_ to make the coupling in the code explicit

mschwoer · 2026-07-13T14:38:50Z

+        m = y_obs.shape[0]
+        df = m - rank
+
+        # DF check: if df <= 0, skip this feature


please check for redundant comments

mschwoer · 2026-07-13T14:40:56Z

+    n_contrasts, P = log2fcs.shape
+
+    # Remove features that have no fit at all, i.e. NaN columns in log2fcs, stdevs_unscaled, or nan in sigma2, or dfs
+    valid = np.isfinite(sigma2) & np.isfinite(dfs) & (dfs > 0)


valid_mask?

mschwoer · 2026-07-13T14:43:45Z

+    )
+    fit = limma.eBayes(fit)
+
+    # scatter back to full (88, P), excluded features stay NaN


(88, P) ?

vbrennsteiner added 4 commits June 26, 2026 17:12

add initial implementation for nan-aware empirical bayes moderated t-…

cc62bcf

…test

output standardization draft

00cfd1b

working, formatted implementation for new diff_exp_ebayes prior to te…

a423f33

…sting implementation#

add comprehensive tests for new differential expression module with a…

646723b

… special focus on preventing any misalignments during the nan_lfit step, wehre the matrices with skipped coefficients are scattered back to full size

vbrennsteiner requested review from lucas-diedrich, mschwoer and shanibmo June 26, 2026 20:52

vbrennsteiner self-assigned this Jun 26, 2026

vbrennsteiner added 4 commits June 28, 2026 22:29

split nan_lfit into one function to make the actual feature mask base…

2ba9104

…d on control samples, and nan_lmfit that purely handles linear fitting on the adata with a design matrix and an optional feature mask

simplify interface by changing high and low filter both to min required

9576b0d

streamline interface by making missingness gates explicitly control/t…

a09334e

…reatment

vbrennsteiner added 9 commits June 29, 2026 12:24

remove pre-gate filtering in favor of symmetric post-fit gating

e0e5811

add functionality forcount-based completeness filtering. Sequester so…

f47b575

…me functionality to small helper functions to make filter_data_completeness less complex

Merge branch 'filter_max_missing_by_count' into expanded_ebayes_ttest

2f8d7d2

add observed=False to groupby operation in filter_data_completeness g…

dbc4c0f

…roupby operation

Merge branch 'filter_max_missing_by_count' into expanded_ebayes_ttest

9af4332

add optional coefficeint return for covariate vs variate diagnostics …

1ec2e24

…and cover with unit tests

add optional coefficeint return for covariate vs variate diagnostics …

a1387cd

…and cover with unit tests

Merge branch 'main' into expanded_ebayes_ttest

f928b48

vbrennsteiner added enhancement New feature or request breaking-change Breaking change that changes the API labels Jul 3, 2026

vbrennsteiner marked this pull request as ready for review July 3, 2026 10:44

mschwoer reviewed Jul 13, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Expanded ebayes ttest#240

Expanded ebayes ttest#240
vbrennsteiner wants to merge 17 commits into
mainfrom
expanded_ebayes_ttest

vbrennsteiner commented Jun 26, 2026 •

edited

Loading

Uh oh!

codecov-commenter commented Jun 29, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

mschwoer Jul 13, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

		return dm, {"condition_col_idxs": condition_col_idxs, "covariate_col_idxs": covariate_col_idxs}


		def summarize_design_matrix(

Uh oh!

Conversation

vbrennsteiner commented Jun 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Purpose of this PR

Critical tests

Outlook

Uh oh!

codecov-commenter commented Jun 29, 2026

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

vbrennsteiner commented Jun 26, 2026 •

edited

Loading