add ont argument to enricher

[YanisChrys]

Hello,

I would like to suggest to add the ont argument from enrichGO to enricher which would make the output immediately compatible with other functions, like simplify. This would help a lot!

Thanks!

[guidohooiveld]

To me that doesn't make that much sense. The whole purpose of the function enricher is that a user can provide a custom 2-column data.frame as input; nothing more, nothing less.

Based on such 2-column data.frame, how would you determine the Gene Ontology category?
If you would like to analyze only a subset of Gene Ontology, I would suggest you do this before 'feeding' the data.frame to enricher(), for example as per the code below.

> ## load libraries
> library(clusterProfiler)
> library(enrichplot)
> library(org.Hs.eg.db)
> library(GO.db)
> 
> 
> ## create TERM2GENE and TERM2NAME for GO
> ## code is based on clusterProfiler code
> ## https://github.com/YuLab-SMU/clusterProfiler/blob/b24be8c372f947ed79f8f42c1a7f6624c14fb125/R/enrichGO.R#L156
> 
> 
> ## set GO category to be analyzed,
> ## i.e. 'BP', 'CC', 'MF', or 'ALL'
> 
> ont <- "BP"  ##limit sets to GO-BP
> 
> goterms <- AnnotationDbi::Ontology(GO.db::GOTERM)
> if (ont != "ALL") {goterms <- goterms[goterms == ont]}
> 
> term2gene.go <- AnnotationDbi::mapIds(org.Hs.eg.db, keys=names(goterms),
+                                       column="ENTREZID", keytype="GOALL", multiVals='list')
'select()' returned 1:many mapping between keys and columns
> term2gene.go <- stack(term2gene.go)
> term2gene.go <- term2gene.go[!is.na(term2gene.go[,"values"]),] [,c(2,1)]
> colnames(term2gene.go) <- c("from","to")
> 
> term2name.go <- Term(GOTERM)[ names(goterms) ]
> term2name.go <- data.frame("from"=names(term2name.go), "to"=term2name.go)
> 
> 
> head(term2gene.go)
        from   to
3 GO:0000012 1161
4 GO:0000012 2074
5 GO:0000012 3981
6 GO:0000012 7014
7 GO:0000012 7141
8 GO:0000012 7374
> head(term2name.go)
                 from                                  to
GO:0000001 GO:0000001           mitochondrion inheritance
GO:0000011 GO:0000011                 vacuole inheritance
GO:0000012 GO:0000012          single strand break repair
GO:0000017 GO:0000017           alpha-glucoside transport
GO:0000018 GO:0000018     regulation of DNA recombination
GO:0000019 GO:0000019 regulation of mitotic recombination
> 

[YanisChrys]

This would "translate" the go terms to ENTREZIDs, what I am interested in is reducing/summarizing the amount of go terms much like clusterProfiler::simplify or rrvgo::reduceSimMatrix would.
For non-model organisms enricher is the one that works best and I simply want to use its output for simplify and it's only compatible if simplify can see the ontology term is not UNKNOWN. The rest of the dataframe is in the same format as the enrichGO output.

[guidohooiveld]

Aha, I think I got your point; in essence you would like to apply the function simplify() to the output of enricher().

Could you please post which non-model organism you are working with, and where the GO annotations are obtained from? This will allow your use case to become concrete.

Also, adding some code, and pasting a set of (example) input ids would be helpful.

[YanisChrys]

Yes, exactly!

I am working on catfish. I can create an orgdb based on Neoarius graeffei using AnnotationForge::makeOrgPackageFromNCBI.
The GO terms are obtained by mapping the proteome to SwissProt and using the GO terms of the UniProt genes. I complement these (genes not found in SwissProt) with GO Terms from eggnog. It's not perfect but it works.

For example what I am doing right now is:

head(catfish_longannotation)
# A tibble: 6 × 5
  TxName         GO         TERM                           ONTOLOGY Gene 
  <chr>          <chr>      <chr>                          <chr>    <chr>
1 anno1.g8792.t1 GO:0001649 osteoblast differentiation     BP       bglap
2 anno1.g8792.t1 GO:0005509 calcium ion binding            MF       bglap
3 anno1.g8792.t1 GO:0005576 extracellular region           CC       bglap
4 anno1.g8792.t1 GO:0008147 structural constituent of bone MF       bglap
5 anno1.g8792.t1 GO:0031214 biomineral tissue development  BP       bglap
6 anno1.g8792.t1 GO:0032571 response to vitamin K          BP       bglap


> head(catfish_genegain)
[1] "gabrg3" "cntn4"  "szt2"   "rnf219" "cks2"   "wdr91" 

catfish_gainBP_enrich <- clusterProfiler::enricher(
  gene=catfish_genegain,
  pvalueCutoff = 0.05,
  pAdjustMethod = "fdr",
  universe = catfish_all_genes,
  minGSSize = 5,
  maxGSSize = 5000,
  qvalueCutoff = 0.2,
  TERM2GENE = catfish_longannotation %>% filter(ONTOLOGY=="BP") %>% dplyr::select(GO,Gene) %>% distinct(),
  TERM2NAME = catfish_longannotation %>% filter(ONTOLOGY=="BP") %>% dplyr::select(GO,TERM) %>% distinct()
)

catfish_gainBP_enrich@ontology <- 'BP'
catfish_gainBP_enrich@keytype <- 'UNIPROT'
pw_df <- pairwise_termsim(x = catfish_gainBP_enrich,showCategory = 500,method = "JC")
catfish_gainBP_enrich2 <- clusterProfiler::simplify(pw_df, cutoff=0.7, by="p.adjust", select_fun=min, measure = "Wang")

[guidohooiveld]

Thanks.

Would you be able to share the file catfish_longannotation? This would allow me to 'play' with it.

[YanisChrys]

Sure thing! I'm sending you the first 6817 genes (100.000 lines) since it has 1 GO per ont per gene
catfish_longannotation.tsv

[guidohooiveld]

Note that you did not say what issue you had with the code you shared in your last post. Yet, I noticed you work with tibbles, etc. It could be that your issue is related to this. You may call me old-fashioned, but I still work with the 'classical' R-objects, and then your code runs fine.

Since you are working with Gene Ontology, I would like to recommend you also add the indirect GO categories to your catfish gene ids. The function buildGOmap() from the library GOSemSim can be used for this.

Below some (classical) R-code to achieve your goal.

> ## load required libraries
> library(clusterProfiler)
> library(GOSemSim)
> library(GO.db)
> library(enrichplot)
> 
> 
> ## Extract all GO *terms* from GO.db
> go.terms <- AnnotationDbi::Ontology(GO.db::GOTERM)
> 
> ## Extract all GO *names* from GO.db
> go.names <- Term(GOTERM)[ names(go.terms) ]
> 
> 
> 
> ## continue with the provided catfish data 
> cf.GO <- read.delim("catfish_longannotation.tsv")
> 
> ## EXTRA STEP
> ## Add all indirect GO annotations to catfish ids.
> ## In other words, add the GO hierarchy.
> ## For gene identifier, either use 'TxName' or 'Gene'
> ## I prefer to use 'TxName'!
> 
> ## Extract all catfish gene ids + direct GO annottaion.
> cf.GO.ids <- cf.GO[, c("GO","TxName")]
> 
> ## check
> dim(cf.GO.ids)
[1] 100028      2
> 
> head(cf.GO.ids)
          GO         TxName
1 GO:0001649 anno1.g8792.t1
2 GO:0005509 anno1.g8792.t1
3 GO:0005576 anno1.g8792.t1
4 GO:0008147 anno1.g8792.t1
5 GO:0031214 anno1.g8792.t1
6 GO:0032571 anno1.g8792.t1
> 
> tail(cf.GO.ids)
               GO           TxName
100023 GO:0006811 anno1.g172752.t1
100024 GO:0006813 anno1.g172752.t1
100025 GO:0006814 anno1.g172752.t1
100026 GO:0017080 anno1.g172752.t1
100027 GO:0051117 anno1.g172752.t1
100028 GO:1903278 anno1.g172752.t1
> 
> 
> ## add indirect GO annotations
> cf.GO.Map <- buildGOmap(cf.GO.ids)
> 
> 
> ## note that many (indirect) ids have been added!
> dim(cf.GO.Map)
[1] 672557      2
> 
> head(cf.GO.Map)
          GO         TxName
1 GO:0001649 anno1.g8792.t1
2 GO:0005509 anno1.g8792.t1
3 GO:0005576 anno1.g8792.t1
4 GO:0008147 anno1.g8792.t1
5 GO:0031214 anno1.g8792.t1
6 GO:0032571 anno1.g8792.t1
> 
> tail(cf.GO.Map)
                GO           TxName
1235303 GO:0140352 anno1.g172752.t1
1235304 GO:1902305 anno1.g172752.t1
1235305 GO:1902307 anno1.g172752.t1
1235306 GO:1903276 anno1.g172752.t1
1235307 GO:1904062 anno1.g172752.t1
1235308 GO:1904064 anno1.g172752.t1
> 
> ## to these indirect annotations, add column GO subclass and name
> cf.GO.Map[,"ONTOLOGY"] <- go.terms[cf.GO.Map[,"GO"] ]
> cf.GO.Map[,"TERM"] <- go.names[cf.GO.Map[,"GO"] ]
> 
> head(cf.GO.Map)
          GO         TxName ONTOLOGY                           TERM
1 GO:0001649 anno1.g8792.t1       BP     osteoblast differentiation
2 GO:0005509 anno1.g8792.t1       MF            calcium ion binding
3 GO:0005576 anno1.g8792.t1       CC           extracellular region
4 GO:0008147 anno1.g8792.t1       MF structural constituent of bone
5 GO:0031214 anno1.g8792.t1       BP  biomineral tissue development
6 GO:0032571 anno1.g8792.t1       BP          response to vitamin K
> 
>

> ## Continue with enricher analysis
> 
> ## As input 250 random genes are used.
> ## No correction for multiple testing is applied (so results are generated!)
> ## Also note that 'TxName' is used for genes,
> ## because 'cf.GO.Map' is 'TxName'-based.
> ## Limit sets to GO-BP.
> 
> ont <- "BP" 
> 
> catfish_gainBP_enrich <- clusterProfiler::enricher(
+   gene = sample( unique( cf.GO[, c("TxName")] ), 250),
+   pvalueCutoff = 0.05,
+   pAdjustMethod = "none",
+   universe = unique( cf.GO[, c("TxName")] ),
+   minGSSize = 5,
+   maxGSSize = 500,
+   qvalueCutoff =1,
+   TERM2GENE = cf.GO.Map[cf.GO.Map[, "ONTOLOGY"] == ont, ] [, c("GO","TxName")],
+   TERM2NAME = cf.GO.Map[cf.GO.Map[, "ONTOLOGY"] == ont, ] [, c("GO","TERM")]
+ )
> 
> 
> ## continue with your code
> catfish_gainBP_enrich@ontology <- ont
> catfish_gainBP_enrich@keytype <- 'UNIPROT'
> pw_df <- pairwise_termsim(x = catfish_gainBP_enrich,
+                           showCategory = 500,
+                           method = "JC")
> 
> catfish_gainBP_enrich2 <- clusterProfiler::simplify(pw_df, cutoff=0.7, by="p.adjust", select_fun=min, measure = "Wang")
> 
> 
> ## check that simplify worked by comparing number of enriched terms; 116 vs 66!
> catfish_gainBP_enrich
#
# over-representation test
#
#...@organism    UNKNOWN 
#...@ontology    BP 
#...@keytype     UNIPROT 
#...@gene        chr [1:250] "anno1.g150403.t1" "anno1.g14792.t1" "anno1.g169303.t1" ...
#...pvalues adjusted by 'none' with cutoff <0.05 
#...116 enriched terms found
'data.frame':   116 obs. of  12 variables:
 $ ID            : chr  "GO:0001822" "GO:0072001" "GO:0061437" "GO:0061440" ...
 $ Description   : chr  "kidney development" "renal system development" "renal system vasculature development" "kidney vasculature development" ...
 $ GeneRatio     : chr  "14/236" "14/236" "4/236" "4/236" ...
 $ BgRatio       : chr  "150/6330" "152/6330" "15/6330" "15/6330" ...
 $ RichFactor    : num  0.0933 0.0921 0.2667 0.2667 0.375 ...
 $ FoldEnrichment: num  2.5 2.47 7.15 7.15 10.06 ...
 $ zScore        : num  3.67 3.61 4.69 4.69 5.04 ...
 $ pvalue        : num  0.00133 0.00151 0.00186 0.00186 0.00249 ...
 $ p.adjust      : num  0.00133 0.00151 0.00186 0.00186 0.00249 ...
 $ qvalue        : num  0.84 0.84 0.84 0.84 0.84 ...
 $ geneID        : chr  "anno1.g14792.t1/anno1.g169303.t1/anno1.g57646.t1/anno1.g56715.t1/anno1.g96989.t1/anno4.16869_t/anno1.g81377.t1/"| __truncated__ "anno1.g14792.t1/anno1.g169303.t1/anno1.g57646.t1/anno1.g56715.t1/anno1.g96989.t1/anno4.16869_t/anno1.g81377.t1/"| __truncated__ "anno1.g169303.t1/anno1.g96989.t1/anno1.g156513.t1/anno1.g105486.t1" "anno1.g169303.t1/anno1.g96989.t1/anno1.g156513.t1/anno1.g105486.t1" ...
 $ Count         : int  14 14 4 4 3 3 3 9 4 4 ...
#...Citation
S Xu, E Hu, Y Cai, Z Xie, X Luo, L Zhan, W Tang, Q Wang, B Liu, R Wang, W Xie, T Wu, L Xie, G Yu. Using clusterProfiler to characterize multiomics data. Nature Protocols. 2024, 19(11):3292-3320 

> 
> catfish_gainBP_enrich2
#
# over-representation test
#
#...@organism    UNKNOWN 
#...@ontology    BP 
#...@keytype     UNIPROT 
#...@gene        chr [1:250] "anno1.g150403.t1" "anno1.g14792.t1" "anno1.g169303.t1" ...
#...pvalues adjusted by 'none' with cutoff <0.05 
#...66 enriched terms found
'data.frame':   66 obs. of  12 variables:
 $ ID            : chr  "GO:0001822" "GO:0042048" "GO:0018107" "GO:0141188" ...
 $ Description   : chr  "kidney development" "olfactory behavior" "peptidyl-threonine phosphorylation" "nucleic acid catabolic process" ...
 $ GeneRatio     : chr  "14/236" "3/236" "3/236" "9/236" ...
 $ BgRatio       : chr  "150/6330" "8/6330" "9/6330" "86/6330" ...
 $ RichFactor    : num  0.0933 0.375 0.3333 0.1047 0.2105 ...
 $ FoldEnrichment: num  2.5 10.06 8.94 2.81 5.65 ...
 $ zScore        : num  3.67 5.04 4.69 3.32 3.99 ...
 $ pvalue        : num  0.00133 0.00249 0.00364 0.0045 0.00469 ...
 $ p.adjust      : num  0.00133 0.00249 0.00364 0.0045 0.00469 ...
 $ qvalue        : num  0.84 0.84 0.84 0.84 0.84 ...
 $ geneID        : chr  "anno1.g14792.t1/anno1.g169303.t1/anno1.g57646.t1/anno1.g56715.t1/anno1.g96989.t1/anno4.16869_t/anno1.g81377.t1/"| __truncated__ "anno2.g58101.t1/anno4.16869_t/anno1.g32197.t1" "anno2.g1210.t1/anno1.g170566.t1/anno2.g4327.t1" "anno1.g184675.t1/anno1.g202228.t1/anno1.g199702.t2/anno1.g81377.t1/anno1.g172770.t1/anno1.g39880.t1/anno1.g8337"| __truncated__ ...
 $ Count         : int  14 3 3 9 4 3 4 10 3 8 ...
#...Citation
S Xu, E Hu, Y Cai, Z Xie, X Luo, L Zhan, W Tang, Q Wang, B Liu, R Wang, W Xie, T Wu, L Xie, G Yu. Using clusterProfiler to characterize multiomics data. Nature Protocols. 2024, 19(11):3292-3320 

> 
> 
> 
> 

[YanisChrys]

Hi Guido,

Thank you for this! It's really useful!

There was nothing technically wrong with my code before (at least I hope) but I just wish that enricher would be compatible with simplify "out of the box" without needing the user to edit the ontology field because the fact that this is needed (and that there is nothing "wrong" with doing that) is not intuitive and initially I thought I couldnt use the two together. But perhaps instead of this requiring a change in enricher maybe it's simplify that needs to have more relaxed input requirements? Not sure.