diff --git a/examples/plot_minimal_pydeseq2_pipeline.py b/examples/plot_minimal_pydeseq2_pipeline.py index 13f1c11d..f2971433 100644 --- a/examples/plot_minimal_pydeseq2_pipeline.py +++ b/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() ` 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 +# `_. + +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: # diff --git a/src/pydeseq2/dds.py b/src/pydeseq2/dds.py index 614bf12d..8bf1d933 100644 --- a/src/pydeseq2/dds.py +++ b/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, ) diff --git a/src/pydeseq2/ds.py b/src/pydeseq2/ds.py index 47c27c84..da7d4dc6 100644 --- a/src/pydeseq2/ds.py +++ b/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 `. """ if coeff not in self.LFC.columns: raise KeyError( diff --git a/tests/test_reference_level.py b/tests/test_reference_level.py new file mode 100644 index 00000000..ba6a43a7 --- /dev/null +++ b/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]"]