<p>Each of these methods do something different, so it is reasonable to expect different results. The bottom line is that there isn't a single question you can ask when you do an &quot;enrichment test&quot;. For example: you can test if a set of genes are more expressed than this other list of genes or they are above the mean levels of any other gene.</p>
<p>Even more, you use the gene ontology, which has an underlying structure (directed acyclic graph) that no method you used takes into account but could modify the results if taken into account. However, the KEGG pathways do not have this structure. So when you want to apply a method mixing both data sources you need to decide if you want to use it or not.</p>
<p>This is also complicated because for most of these methods you can only compare the p-value which depends on the underlying assumptions of the methods (uniform, normal, binomial distribution?)</p>
<p>The latest comprehensive study comparing enrichment methods including some of the methods you used and more is this one (that I know of): <a href="https://doi.org/10.1093/bib/bbz158" rel="nofollow noreferrer">Toward a gold standard for benchmarking gene set enrichment analysis</a>.</p>
<p>Last, what would a correlation tell you here? It is know that there is little agreement between methods. You couldn't say that one method is equal to the other unless you tested with several datasets and you had always a correspondence on the p-value.</p>

CC BY-SA 4.0
2020-03-17T11:34:14.673Z
Thank you, for some reason I was expecting that all the methods would agree on which genesets are differentially enriched to a certain degree.
Darío Rocha
CC BY-SA 4.0
2020-03-17T11:38:19.937Z
Yes, I know it is somehow an expectation of my wet lab colleagues.
llrs
CC BY-SA 4.0
2020-03-16T10:44:10.953Z
2020-11-04T16:19:31.187Z
2020-11-04T16:19:31.187Z
llrs
<p>I have recently ran in R four geneset enrichment analysis in the same database (TCGA:breast cancer) comparing two intrinsic subtypes. The methods I used were:</p>

<ul>
  <li>MIGSA, that imports mGSZ package and combines it with a SEA algorithm. Using RNA-seq TMM normalized counts as input, because the software supports it. All the parameters are as default.</li>
  <li>mGSZ with standard parameters, using RNA-seq TMM normalized log2 cpm. All parameters as default</li>
  <li>FGSEA preranked, using log fold change as ranking metric, obtained from limma deg from RNA-seq TMM normalized counts. All the parameters are as default.</li>
  <li>GSVA + limma: extracted by-sample enrichment score from RNA-seq TMM normalized log2 cpm, then ran a differential expression analysis with geneset enrichment scores as if they were simple gene expression values. I know this is not a formal approach but I wanted to compare. All the parameters are as default.</li>
</ul>

<p>The genesets tested are the full combined list of KEGG, GO:BP, GO:MF and GO:CC. In all cases I trimmed the genesets to genesets of length &lt;500, which leaves a total of 21723 genesets tested.<br>
Sadly, results are confusing. For starters mGSZ for some reason is cropping the output and giving results for only around 9700 of those genesets -I've already sent an email to the package mainteiner asking about that. But that isn't the biggest problem, the most shocking thing is that none of the analyses performedt has a good correlation to the results of any of the other analyses. Here I plotted the -log10(adjusted p-value + 0.0001) of each analysis agaisnt every other analysis, with lines indicating -log10(0.05)=1.3 and -log10(0.01)=2.<br>
<a href="https://i.sstatic.net/VEGsz.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/VEGsz.png" alt="enter image description here"></a>
I was expecting the results to be somehow different, but here I can't see even something remotely close to a correlation between the different results. It is also evident that MIGSA and GSVA+limma output a lot more significatively enriched genesets than the other two methods.<br>
I have gone through the code time and again, and I can't find anything errors. Also I have searched for studies doing this kind of direct comparison between methods, without luck so far. Has anybody experienced this kind of inconsistencies between gene set enrichment analysis results with different tools? What could I possible be doing wrong?.</p>

CC BY-SA 4.0
2020-03-12T23:26:01.260Z
That's a really weird (and interesting) result. Do you know if each package came with its own library of genesets ? Or are they all using the same ?
dc37
CC BY-SA 4.0
2020-03-13T12:23:15.170Z
@Atpoint Although that is a perfectly valid observation, the main question is not about what could be happening and more about whether anyone has ran into similar results when comparing different methods. I edit the main post adding further details.
Darío Rocha
CC BY-SA 4.0
2020-03-13T12:25:09.807Z
@dc37 Yes, the list of genesets was the same for all the tests, they are KEGG, GO:BP, GO:MF and GO:CC.
Darío Rocha
<r><gsea><go-enrichment>
are GSEA and other geneset enrichment analysis supposed to yield extremely different results between them?
# are GSEA and other geneset enrichment analysis supposed to yield extremely different results between them?

(Asked by: username_0 on 2020-03-12)

I have recently ran in R four geneset enrichment analysis in the same database (TCGA:breast cancer) comparing two intrinsic subtypes. The methods I used were:

* MIGSA, that imports mGSZ package and combines it with a SEA algorithm. Using RNA-seq TMM normalized counts as input, because the software supports it. All the parameters are as default.
* mGSZ with standard parameters, using RNA-seq TMM normalized log2 cpm. All parameters as default
* FGSEA preranked, using log fold change as ranking metric, obtained from limma deg from RNA-seq TMM normalized counts. All the parameters are as default.
* GSVA + limma: extracted by-sample enrichment score from RNA-seq TMM normalized log2 cpm, then ran a differential expression analysis with geneset enrichment scores as if they were simple gene expression values. I know this is not a formal approach but I wanted to compare. All the parameters are as default.

The genesets tested are the full combined list of KEGG, GO:BP, GO:MF and GO:CC. In all cases I trimmed the genesets to genesets of length <500, which leaves a total of 21723 genesets tested.  
Sadly, results are confusing. For starters mGSZ for some reason is cropping the output and giving results for only around 9700 of those genesets -I've already sent an email to the package mainteiner asking about that. But that isn't the biggest problem, the most shocking thing is that none of the analyses performedt has a good correlation to the results of any of the other analyses. Here I plotted the -log10(adjusted p-value + 0.0001) of each analysis agaisnt every other analysis, with lines indicating -log10(0.05)=1.3 and -log10(0.01)=2.  
[![enter image description here](https://i.sstatic.net/VEGsz.png)](https://i.sstatic.net/VEGsz.png)
I was expecting the results to be somehow different, but here I can't see even something remotely close to a correlation between the different results. It is also evident that MIGSA and GSVA+limma output a lot more significatively enriched genesets than the other two methods.  
I have gone through the code time and again, and I can't find anything errors. Also I have searched for studies doing this kind of direct comparison between methods, without luck so far. Has anybody experienced this kind of inconsistencies between gene set enrichment analysis results with different tools? What could I possible be doing wrong?.

**Tags:** r, gsea, go-enrichment
> **Comments:**
> **username_1**: That's a really weird (and interesting) result. Do you know if each package came with its own library of genesets ? Or are they all using the same ?  
> **username_0**: @Atpoint Although that is a perfectly valid observation, the main question is not about what could be happening and more about whether anyone has ran into similar results when comparing different methods. I edit the main post adding further details.  
> **username_0**: @dc37 Yes, the list of genesets was the same for all the tests, they are KEGG, GO:BP, GO:MF and GO:CC.  

---

## 1 answers

### Answer by username_2 on 2020-03-16. Score: 2 | Accepted answer

Each of these methods do something different, so it is reasonable to expect different results. The bottom line is that there isn't a single question you can ask when you do an "enrichment test". For example: you can test if a set of genes are more expressed than this other list of genes or they are above the mean levels of any other gene.

Even more, you use the gene ontology, which has an underlying structure (directed acyclic graph) that no method you used takes into account but could modify the results if taken into account. However, the KEGG pathways do not have this structure. So when you want to apply a method mixing both data sources you need to decide if you want to use it or not.

This is also complicated because for most of these methods you can only compare the p-value which depends on the underlying assumptions of the methods (uniform, normal, binomial distribution?)

The latest comprehensive study comparing enrichment methods including some of the methods you used and more is this one (that I know of): [Toward a gold standard for benchmarking gene set enrichment analysis](https://doi.org/10.1093/bib/bbz158).

Last, what would a correlation tell you here? It is know that there is little agreement between methods. You couldn't say that one method is equal to the other unless you tested with several datasets and you had always a correspondence on the p-value.

> **Comments:**
> **username_0**: Thank you, for some reason I was expecting that all the methods would agree on which genesets are differentially enriched to a certain degree.  
> **username_2**: Yes, I know it is somehow an expectation of my wet lab colleagues.  

---

