docs: explain formula-based reference levels for LFC shrinkage

examples/plot_minimal_pydeseq2_pipeline.py

@@ -298,6 +298,46 @@
 
 print(ds.shrunk_LFCs)  # Will be True only if lfc_shrink() was run.
 
+# %%
+# .. _reference_level_ref:
+#
+# Choosing a reference level for shrinkage
+# """"""""""""""""""""""""""""""""""""""""""""""""""
+#
+# A statistical ``contrast`` can compare any two levels, but
+# :meth:`lfc_shrink() <DeseqStats.lfc_shrink>` requires a single fitted coefficient.
+# For example, with three conditions A, B and C and A as the reference, C versus B
+# is a difference of coefficients, not a column to pass to ``lfc_shrink``.
+# Choose the required reference in the design formula before fitting instead.
+# The deprecated ``ref_level`` argument has no effect.
+#
+# Using our two-condition synthetic data, we choose B as the reference and fit a
+# new dataset to obtain the A-versus-B coefficient. ``C()`` marks a variable as
+# categorical; this also matters when category labels are numbers (e.g. batch IDs).
+# ``contr.treatment(base='B')`` selects the reference category. See Formulaic's
+# `treatment coding documentation
+# <https://matthewwardrop.github.io/formulaic/latest/guides/contrasts/#treatment-aka-dummy>`_.
+
+dds_ref_b = DeseqDataSet(
+    counts=counts_df,
+    metadata=metadata,
+    design="~ C(condition, contr.treatment(base='B'))",
+    inference=inference,
+)
+dds_ref_b.deseq2()
+
+# %%
+# Formula-based coding changes the coefficient names. Inspect the fitted columns
+# rather than reusing ``"condition[T.B]"`` from the previous design. We also choose
+# the matching A-versus-B contrast, so the Wald results and shrunk LFCs describe
+# the same comparison. Do not replace the design matrix of an already fitted
+# object without refitting its coefficients.
+
+print(dds_ref_b.varm["LFC"].columns)
+ds_A_vs_B = DeseqStats(dds_ref_b, contrast=["condition", "A", "B"], inference=inference)
+ds_A_vs_B.summary()
+ds_A_vs_B.lfc_shrink(coeff="C(condition, contr.treatment(base='B'))[T.A]")
+
 # %%
 # .. _multifactor_ref:
 #

src/pydeseq2/dds.py

@@ -65,7 +65,8 @@ class DeseqDataSet(ad.AnnData):
         Deprecated.
         Continuous factors are now automatically detected from the design, or cast to categorical using the C() operator in the formula. (default: ``None``).
     ref_level
-        Deprecated.
+        Deprecated and ignored. Set the reference level in ``design`` instead, for
+        example ``design="~ C(condition, contr.treatment(base='B'))"``.
     fit_type
         Either ``"parametric"`` or ``"mean"`` for the type of fitting of dispersions to the mean intensity. ``"parametric"``: fit a dispersion-mean relation via a robust gamma-family GLM. ``"mean"``: use the mean of gene-wise dispersion estimates.
         Will set the fit type for the DEA and the vst transformation.
@@ -225,8 +226,10 @@ def __init__(
 
         if ref_level is not None:
             warnings.warn(
-                "ref_level is deprecated and no longer has any effect. It will be"
-                "removed in a future release.",
+                "ref_level is deprecated and no longer has any effect. It will be "
+                "removed in a future release. Set the reference level in the design "
+                "formula instead, for example "
+                "design=\"~ C(condition, contr.treatment(base='B'))\".",
                 DeprecationWarning,
                 stacklevel=2,
             )

src/pydeseq2/ds.py

@@ -324,7 +324,6 @@ def run_wald_test(self) -> None:
             self.statistics.loc[self.dds.new_all_zeroes_genes] = 0.0
             self.p_values.loc[self.dds.new_all_zeroes_genes] = 1.0
 
-    # TODO update this to reflect the new contrast format
     def lfc_shrink(self, coeff: str, adapt: bool = True) -> None:
         """LFC shrinkage with an apeGLM prior :cite:p:`DeseqStats-zhu2019heavy`.
 
@@ -333,11 +332,23 @@ def lfc_shrink(self, coeff: str, adapt: bool = True) -> None:
         Parameters
         ----------
         coeff
-            The LFC coefficient to shrink.
-            Must be one of the columns of the LFC matrix. (default: ``None``).
+            The name of the fitted LFC coefficient to shrink. Must be a column of
+            ``self.LFC`` (also a column of ``dds.obsm["design_matrix"]``), not a
+            contrast vector or a list of factor levels.
         adapt
             Whether to use the MLE estimates of LFC to adapt the prior.
             If False, the prior scale is set to 1. (``default=True``)
+
+        Notes
+        -----
+        Changing the ``contrast`` passed to :class:`DeseqStats` does not create a
+        new fitted coefficient. If the desired comparison is not a single
+        coefficient, choose its reference level in the design formula before
+        fitting a new :class:`~pydeseq2.dds.DeseqDataSet`, for example
+        ``design="~ C(condition, contr.treatment(base='B'))"``. Inspect
+        ``dds.varm["LFC"].columns`` for the resulting coefficient names and use a
+        matching contrast when creating :class:`DeseqStats`. See the
+        :ref:`reference-level example <reference_level_ref>`.
         """
         if coeff not in self.LFC.columns:
             raise KeyError(

tests/test_reference_level.py

@@ -0,0 +1,49 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from pydeseq2.dds import DeseqDataSet
+
+
+@pytest.mark.parametrize("levels", [("A", "B", "C"), (0, 1, 2)])
+def test_formula_reference_level(levels):
+    samples = [f"sample{i}" for i in range(6)]
+    conditions = np.repeat(levels, 2)
+    metadata = pd.DataFrame({"condition": conditions}, index=samples)
+    counts = pd.DataFrame({"gene1": [10, 12, 20, 22, 30, 32]}, index=samples)
+    reference = levels[1]
+
+    dds = DeseqDataSet(
+        counts=counts,
+        metadata=metadata,
+        design=f"~ C(condition, contr.treatment(base={reference!r}))",
+        n_cpus=1,
+    )
+
+    np.testing.assert_array_equal(
+        dds.obsm["design_matrix"],
+        np.column_stack([np.ones(6), conditions == levels[0], conditions == levels[2]]),
+    )
+    # The third level versus the chosen reference is a single fitted coefficient.
+    np.testing.assert_array_equal(
+        dds.contrast(column="condition", baseline=reference, group_to_compare=levels[2]),
+        [0, 0, 1],
+    )
+
+
+def test_ref_level_warning_points_to_formula():
+    samples = [f"sample{i}" for i in range(4)]
+    counts = pd.DataFrame({"gene1": [10, 12, 20, 22]}, index=samples)
+    metadata = pd.DataFrame({"condition": ["A", "A", "B", "B"]}, index=samples)
+
+    with pytest.warns(DeprecationWarning, match=r"design=.*contr\.treatment"):
+        dds = DeseqDataSet(
+            counts=counts,
+            metadata=metadata,
+            design="~condition",
+            ref_level=["condition", "B"],
+            n_cpus=1,
+        )
+
+    # The deprecated argument still does not change the design.
+    assert dds.obsm["design_matrix"].columns.tolist() == ["Intercept", "condition[T.B]"]