<p>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 &quot;GO enrichment&quot; in many publications).</p>
<p>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).</p>
<p>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 &quot;axon guidance&quot; and  &quot;semaphorin interaction&quot; are significant but &quot;Sema4d in semaphorin signaling&quot; is not, you would stop at &quot;semaphorin signaling&quot; removing the result for &quot;axon guidance&quot;.</p>
<p>This is kind of a strong correction, which is guaranteed to lower the multiple testing correction burden (unless you have no significant genesets).</p>
<p>There are however more refined ways to address this:</p>
<p>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 <strong>elim method</strong> and was described in this <a href="https://academic.oup.com/bioinformatics/article/22/13/1600/193669" rel="nofollow noreferrer">publication</a>.</p>
<p>The same publication discusses another method called <strong>weight</strong> which assigns a weight to a significant category based on the significance of its children.</p>
<p>Another method is implemented in the gProfiler tool and introduced in this 2007 <a href="https://academic.oup.com/nar/article/35/suppl_2/W193/2920757" rel="nofollow noreferrer">publication</a>, which is the <strong>g:SCS method</strong>. 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.</p>
<p>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.</p>
<p>Another relevant, more recent publication is the <a href="https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-017-1571-6" rel="nofollow noreferrer">SetRank</a> 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 <code>fgsea</code>), however, I do not know whether this package has methods you can use or recycle.</p>

CC BY-SA 4.0
2023-01-31T13:28:06.910Z
Great answer! Is there a recommended standard package in R to trim DAGs without that one codes that up yourself?
ATpoint
CC BY-SA 4.0
2023-01-31T15:54:39.783Z
@ATpoint IMHO it depends on the type of test. For Over-Representation Analysis (ORA), which is the one that normally uses GO genesets, packages such as `topGO` (Bioconductor) and `gProfiler2` (CRAN) already implement trimming. For GSEA - what the user was asking about - `SetRank` may be one choice, but I have not used it personally. I am partial to `fgsea` for speed, so I did code my own [pruning functions](https://github.com/langleylab/pruning) some time ago. I believe they are a little crude as they remove entire categories based on the DAG, but don't really address overlaps.
gdagstn
CC BY-SA 4.0
2023-02-01T17:28:47.927Z
Thank you! Both answers provide some very good insight in gene set analysis. Trimming based on DAGs seems very straightforward and easy to do given database such as Reactome provides them in various format already. While searching, I also found a package called [GeneSetCluster](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-020-03784-z) which takes GSA output and correct the results for overlapping genes between gene sets. I haven't tried it yet, but this could be an addition the DAG trimming and SetRank.
Kento
CC BY-SA 4.0
2023-01-30T15:30:42.527Z
2023-01-30T15:30:42.527Z
gdagstn
<p>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.</p>

CC BY-SA 4.0
2023-01-27T17:32:57.377Z
2023-01-27T17:32:57.377Z
ATpoint
<p>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.</p>
<p>I am just wondering if there are ways to mitigate it, or it is not a major concern that I should have?</p>

<pathway><gsea>
Is the potential skewing of FDR and ranking of gene set analysis with gene sets with hierarchical structure a concern?
# Is the potential skewing of FDR and ranking of gene set analysis with gene sets with hierarchical structure a concern?

(Asked by: username_0 on 2023-01-27)

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?

**Tags:** pathway, gsea

---

## 2 answers

### Answer by username_1 on 2023-01-30. Score: 2 | Accepted answer

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](https://academic.oup.com/bioinformatics/article/22/13/1600/193669).

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](https://academic.oup.com/nar/article/35/suppl_2/W193/2920757), 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](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-017-1571-6) 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.

> **Comments:**
> **username_2**: Great answer! Is there a recommended standard package in R to trim DAGs without that one codes that up yourself?  
> **username_1**: @ATpoint IMHO it depends on the type of test. For Over-Representation Analysis (ORA), which is the one that normally uses GO genesets, packages such as `topGO` (Bioconductor) and `gProfiler2` (CRAN) already implement trimming. For GSEA - what the user was asking about - `SetRank` may be one choice, but I have not used it personally. I am partial to `fgsea` for speed, so I did code my own [pruning functions](https://github.com/langleylab/pruning) some time ago. I believe they are a little crude as they remove entire categories based on the DAG, but don't really address overlaps.  
> **username_0**: Thank you! Both answers provide some very good insight in gene set analysis. Trimming based on DAGs seems very straightforward and easy to do given database such as Reactome provides them in various format already. While searching, I also found a package called [GeneSetCluster](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-020-03784-z) which takes GSA output and correct the results for overlapping genes between gene sets. I haven't tried it yet, but this could be an addition the DAG trimming and SetRank.  

---

### Answer by username_2 on 2023-01-27. Score: 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.


---

