I am currently running a simple differential gene expression analysis to identify transcriptomic changes between two condition in ~60 patients. This analysis revealed no significantly differentially expressed genes. logFC are symmetrically distributed around 0 and range from -1.5 to 1.5. The adj.P.val ranges from 0.2 to 1 (vast majority is 1).

Nonethless, I decided to run clusterProfiler::GSEA (using hallmark gene sets from msigdb). For that I ranked the all genes according to the logFC that I got from running the differential gene expression analysis. To my surprise, there were plenty of significantly enriched terms (p.val<1e-16). I know this can happen and is also one of the "strengths" of GSEA, because it can detect changes at pathway-level, even if the changes in the expression level of individual genes are not significant, but I did not expect so many enriched pathways.

So I decided to run an additional test and ran differential gene expression analysis between a completely random subset of individuals (disregarding the two initial conditions). The logFC distribution and adj.P.val distribution looked similar to the "real" differential gene expression analsyis (comparing condition A and B). However, again GSEA showed many significantly enriched pathways. This makes me seriously doubt whether I can trust my GSEA results for "real" differential gene expression analysis.

I would greatly appreciate any insights you have on this topic and also maybe an explanation why this can happen and what I can do to prevent it.

EDIT:

Some more informations to illustrate the issue.

Because basically all adj. P values are one, the volcano plot is sort of meaning less. Instead I am showing the distribution of logFC and unadj.P.value

LogFC:

unadjusted P values:

This is how the top enriched pathway looks like. I used the clusterProfiler::GSEA function.

You are absolutely correct. You have nicely illustrated something that is known to some researchers in the gene set testing area but which is not generally appreciated in the bioinformatics community. The pre-ranked version of GSEA produces over-stated statistical significance. It will indeed produce highly significant pathway results even between randomly chosen groups of samples.

The pre-ranked version of GSEA has never been published or recommended by the GSEA authors themselves, presumably for that reason.

The problem is that pre-ranked GSEA treats genes as independent whereas genes in pathways tend to be positively correlated. Even small positive correlations can produce wildly inflated signficance, with true type I error rates around 40% instead of 5%, as we demonstrated here: https://doi.org/10.1093/nar/gks461. We didn't include pre-ranked GSEA in the simulation in our article, but it would have given the same bad results as the other methods that assume genewise independence.

I have been pointing out the problem with pre-ranked GSEA on this forum for nearly 10 years:

• pre-ranked GSEA within R? + Best DESeq2/limma-voom
• GSEA for RNA-seq analysis
• gene set enrichment
• How to get t-statistic values using edgeR?

There is no assumption that any of the genes are significant in a univariate analysis. The test is one of ranking, where groups of genes with higher than expected ranks (under the null) are considered significant. As you can see in your plot, the hypoxia genes are 'piled up' at the top of the ranked list of genes, when the expectation would be a more uniform distribution across the ranks.