Hi there,

I am trying to run pydeseq2 with one single factor but multiple levels. The factor is called "anatomical_cluster" and the levels are the five values anatomical_cluster can take.

As I want to find differentially expressed genes in one anatomical cluster versus. all others - not doing pairwise comparisons - (doing that for each anatomical cluster), I will use a Wald test and wanted to use something similar to what's in these posts:

• https://support.bioconductor.org/p/86347/
• https://support.bioconductor.org/p/118090/

but they're in R and using deseq2.

In those posts they provide a list in contrast and average the logFC of the ref_level (=last argument in contrast) using listValues. In my case, it would be:

res1 <- results(dds, 
contrast = list( c('EOM'),c('masseter','above_tongue','arches','tongue') ), 
listValues = c(1,-1/4)
)

which corresponds to testing the logFC of (1 * EOM) - (1/4 * masseter + 1/4 * above_tongue + 1/4 * arches + 1/4 * tongue)

I tried applying the same to pydeseq2 and here is my code:

# 1. Build the DESeq2 object with the full design model
dds = DeseqDataSet(
    counts=adata.layers["counts"],
    metadata=adata.obs,
    design="~anatomical_cluster",
    refit_cooks=True,
    inference=inference,
)

# 2. Fit the model and compute dispersions
dds.deseq2()

# 3. Perform the Wald Test with the desired contrast
# Asks: Is the L2FC between 'EOM' and the *average* of the other four groups different from 0?
stat_res_wald = DeseqStats(
    dds, 
    contrast=['anatomical_cluster', 'EOM', ['masseter', 'above_tongue', 'arches', 'tongue']], 
    inference=inference
)

# 4. Print results
stat_res_wald.summary()

This is not working, it says that I'm not providing an existing category of anatomical_cluster. So I wonder: can we provide a list inside 'contrast' argument using pydeseq2? can we adjust for the log2FC like it is done with listValues?

Thank you very much for any precious help!