Is the potential skewing of FDR and ranking of gene set analysis with gene sets with hierarchical structure a concern?

Question

When performing GSEA and many other gene set analysis, when gene sets with hierarchy such as KEGG collection or Reactome collection, the FDR could be skewed sometimes, and the ranking of gene sets by any metrics could be confusing because multiple subsets belonging to the same set would have similar ranking, pushing other gene sets down the list. For example, in reactome axon guidance (R-HSA-422475.6) contains semaphorin interaction (R-HSA-373755.1), which contains Sema4d in semaphorin signaling (R-HSA-400685.2). The enrichment of Sema4d gene set usually (though not necessarily) entails the enrichment of the other gene sets above it. That might lead to an over-adjustment of p-value if I am not mistaken.

I am just wondering if there are ways to mitigate it, or it is not a major concern that I should have?

Answer 1

This is actually a known problem in the correction of nominal p-values in over-representation analyses that use the Gene Ontology categories (commonly but misleadingly called "GO enrichment" in many publications).

Most currently used methods for p-value adjustment (like BH) assume that p-values are independent, which is not true for a case like the one you describe - a set contained by another set is basically by definition not independent. Even if you didn't have a hierarchical structure, you would still have many overlaps across all your genesets, which would still mean that you cannot assume independence. So in both cases, but I think especially in your case, a BH/FWER correction would be statistically wrong (and possibly overly conservative).

There are a few ways to avoid this: one way would be to prune the directed acyclic graph of your genesets/categories, retaining only the categories that are enriched and have no enriched children. In your case, for instance, if "axon guidance" and "semaphorin interaction" are significant but "Sema4d in semaphorin signaling" is not, you would stop at "semaphorin signaling" removing the result for "axon guidance".

This is kind of a strong correction, which is guaranteed to lower the multiple testing correction burden (unless you have no significant genesets).

There are however more refined ways to address this:

One of them is to remove from set A (parent) any genes that belong to a significantly over-represented/enriched set A' (child). This method is termed the elim method and was described in this publication.

The same publication discusses another method called weight which assigns a weight to a significant category based on the significance of its children.

Another method is implemented in the gProfiler tool and introduced in this 2007 publication, which is the g:SCS method. Namely it calculates p-value thresholds based on empirical, randomized p-values for sets of different sizes. According to the authors it is more conservative than BH but less than Bonferroni.

The only issue with these approaches is that there are, to my knowledge, no ways to apply them to your data independently, i.e. they are coded within the respective tools and do not exist as stand-alone packages for multiple test correction. I hope I'm wrong and someone developed these versions, but for the time being the only thing you can do is to try to implement your own versions of these methods.

Another relevant, more recent publication is the SetRank one, which describes a method to do GSEA with a p-value calculation that accounts for overlaps, and has its own package on CRAN. If you prefer to apply a correction to your own results (like the ones you get from fgsea), however, I do not know whether this package has methods you can use or recycle.

Answer 2

Statisticians might kill me for it but I've always filtered long and redundant geneset lists to only contain what remotely makes biological sense in my setup. If I am investigating leukocyte transcriptomes I do not care about terms related to filtration processes in the kidney or signal transduction in the hypothalamus. I think the redundancy of common pathway databases in terms of genes contributing to many terms justifiers to eliminate certain terms and choose those that make sense and are not overly large (for example not more than 500 genes). Top-of-hierarchy terms like 'metabolism' eith more than 1000 genes are anyway unspecific and hard to interpret. That then reduces the multiple testing burden quite a lot.