Issue defining "universe" in enricher function

[GretchenSeim]

Currently trying to perform an enrichment using enricher and a self-curated TERM2GENE list and self defined universe of genes. However even though I am specifying a universe it seems to be pulling the default universe - aka the list of genes in the TERM2GENE list. I have ~8000 genes in my universe and it is just pulling the 30 or so in m TERM2GENE list. Wondering if there is something I'm doing wrong or if there is some sort of bug. There isn't any sort of error message - I can just see in the enrichment result that the universe is much smaller than the 8000 or so genes I'd like it to be.

##defining my universe
universe2<-rownames(DEs)

this results in a character vector of 8563

then I run the enricher as such

kk2<-enricher(gene=rownames(cluster2), pvalueCutoff = 0.05, universe=universe2, TERM2GENE=TERM2GENE)

and my enrichment result shows a universe of only 52 - the number of genes in my TERM2GENE list

[apredeus]

I have similar problem - I think perhaps "universe" definition is working in a way I don't understand?

I've defined universe parameter to all protein-coding genes in my annotation (about 20k symbols), and then ran enrichment vs. MsigDB Hallmark dataset. However, the background was only 4313 genes - exactly the size of the intersection between the "universe" vector and all symbols in H database.

This is changing p-values a lot, so it would be nice to sort it out.

[Jeff-Gui]

Yes I agree with both of you. The universe is just an intersection between the specified vector and all genes in the "TERM2GENE" list. My guess is that the null distribution of p values are generated from randomly sampling the universe genes and then computing the overlap to genes belonging to a term.

[Weihankk]

I also have the same problem. I have ~27K genes as universe genes, ~25k of them have GO annotation terms and as TERM2GENE, then I choose ~10k interesting genes to conduct run enricher. However, the denominator of bgratio is ~25k, and the denominator of generatio is ~ 9k which is the intersect of ~10k interesting genes and ~25k have GO terms. May I ask if this is normal and correct?

Weihan

[avishai987]

Maybe it's more accurate to limit the background to genes that are mapped to the database?
like some say here: https://www.biostars.org/p/439428/

[RaymondSHANG]

Hi, Prof.Yu,
I do have the same issue when using enricher to do the over presentation analysis. The current enricher/DOSE::enricher_internal() used modified hypergeometric test with genes in the (TERM2GENE) as background. This is debatable as discussed:

1. https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1009935
2. https://academic.oup.com/bib/article/22/1/545/5722384

Would be great if we could define our own bg/universe genes in the new versions, rather than the intersect of ourbg with unique(TERM2GENE$gene). If there are really something we need to trim, perhaps we could trim TERM2GENE, removing genes not appeared in the universe, rather than the opposite.
As a clusterProfile User, we occasional analyze pathway level expressions using our customized genesets, which might only contain several pathways, and hundreds of genes. universe defined as all genes in DEG analysis, rather than that in TERM2GENE, make more sense.

Best regards,

[huerqiang]

@RaymondSHANG Please read the doc of enricher():

universe	background genes. If missing, the all genes listed in the database (eg TERM2GENE table) will be used as background.

[RaymondSHANG]

Hi,@huerqiang,
Thanks for the reply. I checked the R functions for enricher().
enricher() --> enricher_internal --> DOSE:::enricher_internal(). And, in DOSE:enricher_internal.R, lines 61-68:

if(!is.null(universe)) {
        if (is.character(universe)) {
            extID <- intersect(extID, universe)
        } else {
            ## https://github.com/YuLab-SMU/clusterProfiler/issues/217
            message("`universe` is not in character and will be ignored...")
        }
}

And line 100:

N <- rep(length(extID), length(M))

When you got intersect(extID, universe), even though universe is provided, only the intersect is selected as the background. Based on my current understanding, N should be:

N <- rep(length(universe), length(M))

I think all above posts in the issue are talking about this.

Thank you,

[RaymondSHANG]

I think recalculating N in DOSE will solve the problem. N <- rep(length(universe), length(M)) TERM2GENE is a kind of global variable, change TERM2GENE may result changes elsewhere in the whole package.

…

On Thu, Oct 20, 2022 at 9:46 AM Graeme Grimes ***@***.***> wrote: Would adding all the genes in the universe as a term in the term2gene database fix this? e.g. # define a TERM2GENE data.frame then add another term including all genes in universe t2uni<-data.frame(term="universe",gene=universe) TERM2GENE=rbind(TERM2GENE,t2uni) — Reply to this email directly, view it on GitHub <#283 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AE6GPCGYR5QDEUMPMXA7QVLWEFZPLANCNFSM4O3IJJ5Q> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- *Raymond "Yuan" SHANG, Ph.D.*

[ariadnaaterrades]

Hi,

We have also faced with this problem. From our point of view, it would be correct to filter the genesets (gmt) to obtain those genes which are present also in your universe. On the other hand, filtering the universe with the genes in the genesets (intersection between universe and gmt) will greatly affect the p value obtained due to the fact that you are reducing the universe (N) to the geneset collection (eg. GO has many annotated genes ~25K, but theoretically a gmt could have only one geneset). A quick solution for the time being, not the best one though, could be adding your universe to your gmt as a pathway, so all the genes in your universe will be part of the term2gene and the intersection will not affect your universe size nor your significance level.

It would be nice to find a solution to be able to avoid obtaining underestimated significance.

Thank you,
Ariadna

[christianmaueroeder]

Hi all,

I am experiencing the same problem at the moment. Is this going to be fixed or is it necessary to adjust manually, as suggested by @RaymondSHANG ?

How did you proceed in this case?

Cheers

[SamGG]

@RaymondSHANG I led to the same conclusion. Do you plan a PR to DOSE?

[GuangchuangYu]

If you want to keep the universe untouched, aka, not performing intersection with genes that have annotations, you can force it to do so via options(enrichment_force_universe = TRUE).

Please try it and give me feedback if there is something not work as expected.

[christianmaueroeder]

If you want to keep the universe untouched, aka, not performing intersection with genes that have annotations, you can force it to do so via options(enrichment_force_universe = TRUE).

Please try it and give me feedback if there is something not work as expected.

Thanks for this option, however, I think what is done now is taking all genes that are stored in OrgDb irrespective of the values specified by universe. For example, after using org.Mm.eg.db, I get a background of 28891 genes, thus, all mouse genes and not only the genes that are expressed in my target cells.