5
$\begingroup$

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 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?.

$\endgroup$
4
  • $\begingroup$ 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 ? $\endgroup$ – dc37 Mar 12 '20 at 23:26
  • 3
    $\begingroup$ Without code and details on how exactly you ran the tools (and with which input) this is hard to answer. $\endgroup$ – ATpoint Mar 12 '20 at 23:26
  • $\begingroup$ @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. $\endgroup$ – Darío Rocha Mar 13 '20 at 12:23
  • $\begingroup$ @dc37 Yes, the list of genesets was the same for all the tests, they are KEGG, GO:BP, GO:MF and GO:CC. $\endgroup$ – Darío Rocha Mar 13 '20 at 12:25
2
$\begingroup$

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.

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.

$\endgroup$
2
  • $\begingroup$ Thank you, for some reason I was expecting that all the methods would agree on which genesets are differentially enriched to a certain degree. $\endgroup$ – Darío Rocha Mar 17 '20 at 11:34
  • $\begingroup$ Yes, I know it is somehow an expectation of my wet lab colleagues. $\endgroup$ – llrs Mar 17 '20 at 11:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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