I am working with an RNA-seq dataset generated from a complex experimental design. Our main question is focused on the interaction between two variables, one (INJ) with two conditions (L,S) and the other (SOC) with four conditions (A,B,C,D). We also have a third variable for tissue type (TIS) from 3 tissues (H,N,T), although we are not interested in this variable for this particular question.

We were able extract DE genes underlying the interaction of interest in DESeq2 using an LRT test. We also ran WGCNA to get gene clusters and looked at trait correlations. Using these WGNCA clusters and VST normalized count data, I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables.

What I've tried so far was generating a pseudo-variable that is a combination of INJ and SOC with 8 separate conditions for running in GSEA, but this doesn't quite get at the interaction between variables since we really want to see if the effect of INJ differs between SOC conditions.

One thing I was thinking was running a GSEA-like analysis on log2 ratios between INJ states for each SOC condition, but these would have to be summary statistics for groups of samples and I would lose the variance of expression within groups (maybe that doesn't matter though for exploratory analysis?). This would also reduce my "sample" count to 1 per condition, which isn't ideal for this type of analysis and I would probably have to look at differences between gene sets instead of phenotypes for more robust stats.

Any thoughts on ways to do this? I don't need to use GSEA if there is another analysis I could run on gene sets that allows for more complex experimental models. Thanks in advance.

  • $\begingroup$ You could read about differential co-expression. But that would be feasible depending on the number of samples you have for each condition... $\endgroup$
    – llrs
    Commented Sep 29, 2020 at 8:29

1 Answer 1


I'm not sure I understand "I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables". I would run an association analysis (regression) for module eigengenes and the appropriate interaction term. The (most) significant modules are your candidates. This step is simply regression (in the simplest case a linear model should be sufficient), no GSEA needed. You can then subject the identified modules to GSEA or similar to get an idea of their biological function.

  • $\begingroup$ Thanks for the input and sorry for the late response. I've been looking into association analyses for modules/eigengenes and I'm thinking of running LRT analyses using GLMs with and without the interaction term of interest. Would you suggest I run separate analyses using each module's eigengene values individually (vs phenotypes) then correcting for multiple hypothesis, or can you recommend a more elegant way to test all modules at once? We would also like to run 3 separate analyses for each of 20 modules based on tissue sample subsets, so running 60 different analyses is not ideal. Thanks. $\endgroup$
    – jfaberha
    Commented Oct 9, 2020 at 1:39
  • 1
    $\begingroup$ For module eigengenes, a linear model such as ones implemented in limma should be appropriate, no reason to run GLMs. With a linear model, you can simply test the significance of the interaction coefficient which is output by all linear model functions that I know (as long as the full and reduced model only differ by one coefficient). As far as I am aware, the only way to run such association analysis is fitting the models to the eigengenes one by one (perhaps with some empirical Bayes modifications as in limma) and then correcting for multiple comparisons. $\endgroup$ Commented Oct 10, 2020 at 3:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.