Switch to formulaic-contrasts

[grst]

PR Checklist

• Referenced issue is linked
• If you've fixed a bug or added code that should be tested, add tests!
• Documentation in docs is updated

Description of changes

Technical details

Additional context

Close #610

[Zethson]

@grst I'll fix the unrelated tests soon. Sorry about them + the currently uber slow CI.

[grst]

What's left is to update the differential expression tutorial to use a numeric contrast vector instead of a tuple of strings, i.e.

-res_df = pds2.test_contrasts(["Treatment", "Chemo", "Anti-PD-L1+Chemo"])
+res_df = pds2.test_contrasts(pds2.contrast(column="Treatment", baseline="Chemo", group_to_compare="Anti-PD-L1+Chemo"))

I tried a bit, but rerunning the tutorial with edgeR and rpy2 is a huge rabbit hole.

[emdann]

LGTM. I've pushed some edits to pass the rpy2/edgeR tests #692

The main thing to work on seems to be the documentation, e.g. is this still accurate?

pertpy/pertpy/tools/_differential_gene_expression/_base.py

Lines 982 to 992 in 298db0f

	def test_contrasts(self, contrasts, **kwargs):
	"""
	Perform a comparison as specified in a contrast vector.
	Args:
	contrasts: Either a numeric contrast vector, or a dictionary of numeric contrast vectors.
	**kwargs: passed to the respective implementation.
	Returns:
	A dataframe with the results.
	"""

Also typing in base functions like test_contrasts and compare_groups would be helpful.

[emdann]

I'm gonna have a go at editing the tutorial.

[emdann]

One small note: I noticed that the model doesn't complain if you try to specify a complex interaction contrast on a model that wasn't fit with the interaction in the design, but just throws nonsense results.

Following example from the tutorial:

# Exclude patient with progressive disease, or not full rank for interaction
pdata2 = pdata[pdata.obs['Efficacy'] != 'PD'].copy()

# Bad design definition without interaction
pds2 = pt.tl.PyDESeq2(adata=pdata2, design="~ Efficacy + Treatment")
pds2.fit()

interaction_contrast = (
    pds2.cond(Treatment="Chemo", Efficacy="PR") - pds2.cond(Treatment="Chemo", Efficacy="SD")
) - (
    pds2.cond(Treatment="Anti-PD-L1+Chemo", Efficacy="PR") - pds2.cond(Treatment="Anti-PD-L1+Chemo", Efficacy="SD")
)
res_df = pds2.test_contrasts(contrasts=interaction_contrast)

No complaint, but the results are broken:

Log2 fold change & Wald test p-value, contrast vector: [0. 0. 0.]
           baseMean  log2FoldChange  lfcSE  stat  pvalue  padj
A1BG      16.408605             0.0    0.0   NaN     NaN   NaN
A1BG-AS1   1.958737             0.0    0.0   NaN     NaN   NaN
A1CF       0.002053             0.0    0.0   NaN     NaN   NaN
A2M       30.296881             0.0    0.0   NaN     NaN   NaN
A2M-AS1    0.557092             0.0    0.0   NaN     NaN   NaN
...             ...             ...    ...   ...     ...   ...
ZXDC       6.114098             0.0    0.0   NaN     NaN   NaN
ZYG11A     0.093600             0.0    0.0   NaN     NaN   NaN
ZYG11B     3.404941             0.0    0.0   NaN     NaN   NaN
ZYX       77.175203             0.0    0.0   NaN     NaN   NaN
ZZEF1      9.752162             0.0    0.0   NaN     NaN   NaN

While results make sense if I specify the design properly pds2 = pt.tl.PyDESeq2(adata=pdata2, design="~ Efficacy + Treatment + Efficacy*Treatment").

Do we want this to happen? Can we add an informative error if the contrast vector is all zeros?

[Zethson]

<3 looks great!

[Zethson]

@emdann I added type hints as you requested above and updated the submodule to have your tutorial changes as well. With passing CI (I work on it), I am happy to merge this now.

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 86.84211% with 5 lines in your changes missing coverage. Please review.

Project coverage is 64.85%. Comparing base (9bba130) to head (a312bd0).
Report is 5 commits behind head on main.

Files with missing lines	Patch %	Lines
...ertpy/tools/_differential_gene_expression/_base.py	80.00%	3 Missing ⚠️
...y/tools/_differential_gene_expression/_pydeseq2.py	75.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #682      +/-   ##
==========================================
- Coverage   65.56%   64.85%   -0.72%     
==========================================
  Files          47       46       -1     
  Lines        6105     5992     -113     
==========================================
- Hits         4003     3886     -117     
- Misses       2102     2106       +4     

Files with missing lines	Coverage Δ
...py/tools/_differential_gene_expression/__init__.py	100.00% <100.00%> (ø)
...rtpy/tools/_differential_gene_expression/_edger.py	84.84% <100.00%> (-2.25%)	⬇️
...ools/_differential_gene_expression/_statsmodels.py	100.00% <ø> (ø)
pertpy/tools/_distances/_distances.py	89.96% <100.00%> (-0.07%)	⬇️
pertpy/tools/_milo.py	61.75% <100.00%> (ø)
pertpy/tools/_mixscape.py	79.12% <100.00%> (ø)
...y/tools/_differential_gene_expression/_pydeseq2.py	92.30% <75.00%> (-1.03%)	⬇️
...ertpy/tools/_differential_gene_expression/_base.py	24.84% <80.00%> (-7.49%)	⬇️

[Zethson]

Merging now. Making a release as well. @emdann I'll create an issue for your concern above.