pertpy icon indicating copy to clipboard operation
pertpy copied to clipboard

Model doesn't complain if a complex interaction contrast is specified that wasn't fit with the interaction in the design

Open Zethson opened this issue 1 year ago • 0 comments

Report

Raised by @emdann in https://github.com/scverse/pertpy/pull/682#issuecomment-2569982972

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?

Version information

No response

Zethson avatar Jan 04 '25 15:01 Zethson