Gene list functional enrichment analysis and namespace conversion with gprofiler2

Enrichment analysis

gost enables to perform functional profiling of gene lists. The function performs statistical enrichment analysis to find over-representation of functions from Gene Ontology, biological pathways like KEGG and Reactome, human disease annotations, etc. This is done with the hypergeometric test followed by correction for multiple testing.

A standard input of the gost function is a (named) list of gene identifiers. The list can consist of mixed types of identifiers (proteins, transcripts, microarray IDs, etc), SNP IDs, chromosomal intervals or functional term IDs.

The parameter organism enables to define the corresponding source organism for the gene list. The organism names are usually constructed by concatenating the first letter of the name and the family name, e.g human - hsapiens. If some of the input gene identifiers are fully numeric, the parameter numeric_ns enables to define the corresponding namespace. See section Supported organisms and identifier namespaces for links to supported organisms and namespaces.

If the input genes are decreasingly ordered based on some biological importance, then ordered_query = TRUE will take this into account. For instance, the genes can be ordered according to differential expression or absolute expression values. In this case, incremental enrichment testing is performed with increasingly larger numbers of genes starting from the top of the list. Note that with this parameter, the query size might be different for every functional term.

The parameter significant = TRUE is an indicator whether all or only statistically significant results should be returned.

In case of Gene Ontology (GO), the exclude_iea = TRUE would exclude the electronic GO annotations from the data source before testing. These are the terms with the IEA evidence code indicating that these annotations are assigned to genes using in silico curation methods.

In order to measure under-representation instead of over-representation set measure_underrepresentation = TRUE.

By default, the user_threshold = 0.05 which defines a custom p-value significance threshold for the results. Results with smaller p-value are tagged as significant. We don’t recommend to set it higher than 0.05.

In order to reduce the amount of false positives, a multiple testing correction method is applied to the enrichment p-values. By default, our tailor-made algorithm g:SCS is used (correction_method = "gSCS" with synonyms g_SCS and analytical), but there are also options to apply the Bonferroni correction (correction_method = "bonferroni") or FDR (correction_method = "fdr"). The adjusted p-values are reported in the results.

The parameter domain_scope defines how the statistical domain size is calculated. This is one of the parameters in the hypergeometric probability function. If domain_scope = "annotated" then only the genes with at least one annotation are considered to be part of the full domain. In case if domain_scope = "known" then all the genes of the given organism are considered to be part of the domain.

Depending on the research question, in some occasions it is advisable to limit the domain/background set. For example, one may use the custom background when they want to compare a gene list with a custom list of expressed genes. gost provides the means to define a custom background as a (mixed) vector of gene identifiers with the parameter custom_bg. If this parameter is used, then the domain scope is set to domain_scope = "custom". It is also possible to set this parameter to domain_scope = "custom_annotated" which will use the set of genes that are annotated in the data source and are also included in the user provided background list.

The parameter sources enables to choose the data sources of interest. By default, all the sources are analysed. The available data sources and their abbreviations are listed under section [Data sources]. For example, if sources = c("GO:MF", "REAC") then only the results from molecular functions branch of Gene Ontology and the pathways from Reactome are returned. One can also upload their own annotation data which is further described in the section Custom data sources with upload_GMT_file.

Parameter highlight includes an indicator TRUE/FALSE column called “highlighted” to the analysis results to highlight driver terms in GO. This option works starting from Ensembl version 108 (e108). The option doesn’t work with custom GMT files.

gostres <- gost(query = c("X:1000:1000000", "rs17396340", "GO:0005005", "ENSG00000156103", "NLRP1"), 
                organism = "hsapiens", ordered_query = FALSE, 
                multi_query = FALSE, significant = TRUE, exclude_iea = FALSE, 
                measure_underrepresentation = FALSE, evcodes = FALSE, 
                user_threshold = 0.05, correction_method = "g_SCS", 
                domain_scope = "annotated", custom_bg = NULL, 
                numeric_ns = "", sources = NULL, as_short_link = FALSE, highlight = TRUE)

The result is a named list where “result” is a data.frame with the enrichment analysis results and “meta” containing a named list with all the metadata for the query.

names(gostres)
#> [1] "result" "meta"

head(gostres$result, 3)
#>     query significant      p_value term_size query_size intersection_size
#> 1 query_1        TRUE 4.998799e-02         1          3                 1
#> 2 query_1        TRUE 2.115706e-30        51         20                14
#> 3 query_1        TRUE 7.617347e-18       249         20                13
#>   precision     recall    term_id source                         term_name
#> 1 0.3333333 1.00000000 CORUM:6180  CORUM           PPP2R1A-PPP2R3B complex
#> 2 0.7000000 0.27450980 GO:0048013  GO:BP ephrin receptor signaling pathway
#> 3 0.6500000 0.05220884 GO:0097485  GO:BP        neuron projection guidance
#>   effective_domain_size source_order                parents highlighted
#> 1                  3383         1952          CORUM:0000000       FALSE
#> 2                 21031        12892             GO:0007169        TRUE
#> 3                 21031        19215 GO:0031175, GO:0048812       FALSE

The result data.frame contains the following columns:

• query - the name of the input query which by default is the order of query with the prefix “query_”. This can be changed by using a named list input.
• significant - indicator for statistically significant results
• p_value - hypergeometric p-value after correction for multiple testing
• term_size - number of genes that are annotated to the term
• query_size - number of genes that were included in the query. This might be different from the size of the original list if:
  1. any genes were mapped to multiple Ensembl gene IDs
  2. any genes failed to be mapped to Ensembl gene IDs
  3. the parameter ordered_query = TRUE and the optimal cutoff for the term was found before the end of the query
  4. the domain_scope was set to “annotated” or “custom”
• intersection_size - the number of genes in the input query that are annotated to the corresponding term
• precision - the proportion of genes in the input list that are annotated to the function (defined as intersection_size/query_size)
• recall - the proportion of functionally annotated genes that the query recovers (defined as intersection_size/term_size)
• term_id - unique term identifier (e.g GO:0005005)
• source - the abbreviation of the data source for the term (e.g. GO:BP)
• term_name - the short name of the function
• effective_domain_size - the total number of genes “in the universe” used for the hypergeometric test
• source_order - numeric order for the term within its data source (this is important for drawing the results)
• parents - list of term IDs that are hierarchically directly above the term. For non-hierarchical data sources this points to an artificial root node
• highlighted - TRUE/FALSE for GO terms to indicate driver terms detected by a two-stage algorithm for filtering (works starting from version e108). NB! For non-GO terms the value is always FALSE. This column could be used to indicate terms to highlight in the gostplot (see below).

names(gostres$meta)
#> [1] "query_metadata"  "result_metadata" "genes_metadata"  "timestamp"      
#> [5] "version"

The query parameters are listed in the “query_metadata” part of the metadata object. The “result_metadata” includes the statistics of data sources that are used in the enrichment testing. This includes the “domain_size” showing the number of genes annotated to this domain. The “number_of_terms” indicating the number of terms g:Profiler has in the database for this source and the nominal significance “threshold” for this source. The “genes_metadata” shows the specifics of the query genes (failed, ambiguous or duplicate inputs) and their mappings to the ENSG namespace. In addition, the query time and the used g:Profiler data version are shown in the metadata.

The parameter evcodes = TRUE includes the evidence codes to the results. In addition, a column “intersection” will appear to the results showing the input gene IDs that intersect with the corresponding functional term. Note that his parameter can decrease the performance and make the query slower.

gostres2 <- gost(query = c("X:1000:1000000", "rs17396340", "GO:0005005", "ENSG00000156103", "NLRP1"), 
                organism = "hsapiens", ordered_query = FALSE, 
                multi_query = FALSE, significant = TRUE, exclude_iea = FALSE, 
                measure_underrepresentation = FALSE, evcodes = TRUE, 
                user_threshold = 0.05, correction_method = "g_SCS", 
                domain_scope = "annotated", custom_bg = NULL, 
                numeric_ns = "", sources = NULL, highlight = TRUE)

head(gostres2$result, 3)
#>     query significant      p_value term_size query_size intersection_size
#> 1 query_1        TRUE 4.998799e-02         1          3                 1
#> 2 query_1        TRUE 2.115706e-30        51         20                14
#> 3 query_1        TRUE 7.617347e-18       249         20                13
#>   precision     recall    term_id source                         term_name
#> 1 0.3333333 1.00000000 CORUM:6180  CORUM           PPP2R1A-PPP2R3B complex
#> 2 0.7000000 0.27450980 GO:0048013  GO:BP ephrin receptor signaling pathway
#> 3 0.6500000 0.05220884 GO:0097485  GO:BP        neuron projection guidance
#>   effective_domain_size source_order                parents
#> 1                  3383         1952          CORUM:0000000
#> 2                 21031        12892             GO:0007169
#> 3                 21031        19215 GO:0031175, GO:0048812
#>                                                                                                                                    evidence_codes
#> 1                                                                                                                                           CORUM
#> 2 IDA IBA IEA,ISS IBA IEA,IBA IEA,IDA IBA,IGI IBA IEA,ISS IBA IEA,IDA IBA IEA,IGI ISS IBA IEA,ISS IBA IEA,IBA,IDA IBA IEA,IDA IBA,IBA IEA,IBA IEA
#> 3                                         IBA,ISS IBA IEA,IBA,IBA IEA,ISS IBA IEA,ISS IBA IEA,IBA IEA,IBA,ISS IBA,ISS IBA IEA,ISS IBA IEA,IBA,IBA
#>                                                                                                                                                                                                                      intersection
#> 1                                                                                                                                                                                                                 ENSG00000167393
#> 2 ENSG00000044524,ENSG00000070886,ENSG00000080224,ENSG00000108947,ENSG00000116106,ENSG00000133216,ENSG00000135333,ENSG00000143590,ENSG00000145242,ENSG00000146904,ENSG00000154928,ENSG00000182580,ENSG00000183317,ENSG00000243364
#> 3                 ENSG00000044524,ENSG00000070886,ENSG00000080224,ENSG00000108947,ENSG00000116106,ENSG00000133216,ENSG00000135333,ENSG00000143590,ENSG00000145242,ENSG00000154928,ENSG00000182580,ENSG00000183317,ENSG00000243364
#>   highlighted
#> 1       FALSE
#> 2        TRUE
#> 3       FALSE

The result data.frame will include additional columns:

• evidence_codes - a lists of all evidence codes for the intersecting genes between input and the term. The evidences are separated by comma for each gene.
  
• intersection - a comma separated list of genes from the query that are annotated to the corresponding term

The query results can also be gathered into a short-link to the g:Profiler web tool. For that, set the parameter as_short_link = TRUE. In this case, the function gost() returns only the web tool link to the results as a character string. For example, this is useful when you discover an interesting result you want to instantly share with your colleagues. Then you can just programmatically generate the short-link and copy it to your colleagues.

gostres_link <- gost(query = c("X:1000:1000000", "rs17396340", "GO:0005005", "ENSG00000156103", "NLRP1"), 
                as_short_link = TRUE)

This query returns a short-link of form https://biit.cs.ut.ee/gplink/l/HfapQyB5TJ.

Multiple queries

The function gost also allows to perform enrichment on multiple input gene lists. Multiple queries are automatically detected if the input query is a list of vectors with gene identifiers and the results are combined into identical data.frame as in case of single query.

multi_gostres1 <- gost(query = list("chromX" = c("X:1000:1000000", "rs17396340", 
                                                 "GO:0005005", "ENSG00000156103", "NLRP1"),
                             "chromY" = c("Y:1:10000000", "rs17396340", 
                                          "GO:0005005", "ENSG00000156103", "NLRP1")), 
                       multi_query = FALSE)

head(multi_gostres1$result, 3)
#>    query significant      p_value term_size query_size intersection_size
#> 1 chromX        TRUE 4.998799e-02         1          3                 1
#> 2 chromX        TRUE 2.115706e-30        51         20                14
#> 3 chromX        TRUE 7.617347e-18       249         20                13
#>   precision     recall    term_id source                         term_name
#> 1 0.3333333 1.00000000 CORUM:6180  CORUM           PPP2R1A-PPP2R3B complex
#> 2 0.7000000 0.27450980 GO:0048013  GO:BP ephrin receptor signaling pathway
#> 3 0.6500000 0.05220884 GO:0007411  GO:BP                     axon guidance
#>   effective_domain_size source_order                parents
#> 1                  3383         1952          CORUM:0000000
#> 2                 21031        12892             GO:0007169
#> 3                 21031         2940 GO:0007409, GO:0097485

The column “query” in the result data.frame will now contain the corresponding name for the query. If no name is specified, then the query name is defined as the order of query with the prefix “query_”.

Another option for multiple gene lists is setting the parameter multiquery = TRUE. Then the results from all of the input queries are grouped according to term IDs for better comparison.

multi_gostres2 <- gost(query = list("chromX" = c("X:1000:1000000", "rs17396340",
                                                 "GO:0005005", "ENSG00000156103", "NLRP1"),
                             "chromY" = c("Y:1:10000000", "rs17396340", 
                                          "GO:0005005", "ENSG00000156103", "NLRP1")), 
                       multi_query = TRUE)

head(multi_gostres2$result, 3)
#>      term_id                   p_values significant term_size query_sizes
#> 1 GO:0005005 1.696883e-42, 2.059426e-35  TRUE, TRUE        14      21, 50
#> 2 GO:0005003 5.185748e-39, 6.260011e-32  TRUE, TRUE        18      21, 50
#> 3 GO:0048013 2.115706e-30, 1.357066e-24  TRUE, TRUE        51      20, 40
#>   intersection_sizes source                              term_name
#> 1             14, 14  GO:MF transmembrane-ephrin receptor activity
#> 2             14, 14  GO:MF               ephrin receptor activity
#> 3             14, 14  GO:BP      ephrin receptor signaling pathway
#>   effective_domain_size source_order    parents
#> 1                 20212         1236 GO:0005003
#> 2                 20212         1234 GO:0004714
#> 3                 21031        12892 GO:0007169

The result data.frame contains the following columns:

• term_id - unique term identifier (e.g GO:0005005)
• p_values - hypergeometric p-values after correction for multiple testing in the order of input queries
• significant - indicators in the order of input queries for statistically significant results
• term_size - number of genes that are annotated to the term
• query_sizes - number of genes that were included in the query in the order of input queries
• intersection_sizes - the number of genes in the input query that are annotated to the corresponding term in the order of input queries
• source - the abbreviation of the data source for the term (e.g. GO:BP)
• term_name - the short name of the function
• effective_domain_size - the total number of genes “in the universe” used for the hypergeometric test
• source_order - numeric order for the term within its data source (this is important for drawing the results)
• source_order - numeric order for the term within its data source (this is important for drawing the results)
• parents - list of term IDs that are hierarchically directly above the term. For non-hierarchical data sources this points to an artificial root node.

Visualization

A major update in this package is providing the functionality to produce similar visualizations as are now available from the web tool.

The enrichment results are visualized with a Manhattan-like-plot using the function gostplot and the previously found gost results gostres:

gostplot(gostres, capped = TRUE, interactive = TRUE)

The x-axis represents the functional terms that are grouped and color-coded according to data sources and positioned according to the fixed “source_order”. The order is defined in a way that terms that are closer to each other in the source hierarchy are also next to each other in the Manhattan plot. The source colors are adjustable with the parameter pal that defines the color map with a named list.

The y-axis shows the adjusted p-values in negative log10 scale. Every circle is one term and is sized according to the term size, i.e larger terms have larger circles. If interactive = TRUE, then an interactive plot is returned using the plotly package. Hovering over the circle will show the corresponding information.

The parameter capped = TRUE is an indicator whether the -log10(p-values) would be capped at 16 if bigger than 16. This fixes the scale of y-axis to keep Manhattan plots from different queries comparable and is also intuitive as, statistically, p-values smaller than that can all be summarised as highly significant.

If interactive = FALSE, then the function returns a static ggplot object.

p <- gostplot(gostres, capped = FALSE, interactive = FALSE)
p

The function publish_gostplot takes the static plot object as an input and enables to highlight a selection of interesting terms from the results with numbers and table of results. These can be set with parameter highlight_terms listing the term IDs in a vector or as a data.frame with column “term_id” such as a subset of the result data.frame.

pp <- publish_gostplot(p, highlight_terms = c("GO:0048013", "REAC:R-HSA-3928663"), 
                       width = NA, height = NA, filename = NULL )

The function also allows to save the result into an image file into PNG, PDF, JPEG, TIFF or BMP with parameter filename. The plot width and height can be adjusted with corresponding parameters width and height.

If additional graphical parameters like increased resolution are needed, then the plot objects can easily be saved using the function ggsave from the package ggplot2.

The gost results can also be visualized with a table. The publish_gosttable will create a nice-looking table with the result statistics for the highlight_terms from the result data.frame. The highlight_terms can be a vector of term IDs or a subset of the results.

publish_gosttable(gostres, highlight_terms = gostres$result[c(1:2,10,120),],
                        use_colors = TRUE, 
                        show_columns = c("source", "term_name", "term_size", "intersection_size"),
                        filename = NULL)
#> The input 'highlight_terms' is a data.frame. The column 'term_id' will be used.

The parameter use_colors = FALSE indicates that the p-values column should not be highlighted with background colors. The show_columns is used to list the names of additional columns to show in the table in addition to the “term_id” and “p_value”.

The same functions work also in case of multiquery results showing multiple Manhattan plots on top of each other:

gostplot(multi_gostres2, capped = TRUE, interactive = TRUE)

Note that if a term is clicked on one of the Manhattan plots, it is also highlighted in the others (if it is present) enabling to compare the multiple queries. The insignificant terms are shown with lighter color.

publish_gosttable(multi_gostres1, 
                         highlight_terms = multi_gostres1$result[c(1, 82, 176),],
                        use_colors = TRUE, 
                        show_columns = c("source", "term_name", "term_size"),
                        filename = NULL)
#> The input 'highlight_terms' is a data.frame. The column 'term_id' will be used.

Data sources and versions

Available data sources and their abbreviations are:

• Gene Ontology (GO or by branch GO:MF, GO:BP, GO:CC)
• KEGG (KEGG)
• Reactome (REAC)
• WikiPathways (WP)
• TRANSFAC (TF)
• miRTarBase (MIRNA)
• Human Protein Atlas (HPA)
• CORUM (CORUM)
• Human phenotype ontology (HP)

The function get_version_info enables to obtain the full metadata about the versions of different data sources for a given organism.

get_version_info(organism = "hsapiens")
#> $biomart
#> [1] "Ensembl"
#> 
#> $biomart_version
#> [1] "111"
#> 
#> $display_name
#> [1] "Human"
#> 
#> $genebuild
#> [1] "GRCh38.p14"
#> 
#> $gprofiler_version
#> [1] "e111_eg58_p18_f463989d"
#> 
#> $organism
#> [1] "hsapiens"
#> 
#> $sources
#> $sources$CORUM
#> $sources$CORUM$name
#> [1] "CORUM protein complexes"
#> 
#> $sources$CORUM$version
#> [1] "28.11.2022 Corum 4.1"
#> 
#> 
#> $sources$`GO:BP`
#> $sources$`GO:BP`$name
#> [1] "biological process"
#> 
#> $sources$`GO:BP`$version
#> [1] "annotations: BioMart\nclasses: releases/2024-01-17"
#> 
#> 
#> $sources$`GO:CC`
#> $sources$`GO:CC`$name
#> [1] "cellular component"
#> 
#> $sources$`GO:CC`$version
#> [1] "annotations: BioMart\nclasses: releases/2024-01-17"
#> 
#> 
#> $sources$`GO:MF`
#> $sources$`GO:MF`$name
#> [1] "molecular function"
#> 
#> $sources$`GO:MF`$version
#> [1] "annotations: BioMart\nclasses: releases/2024-01-17"
#> 
#> 
#> $sources$HP
#> $sources$HP$name
#> [1] "Human Phenotype Ontology"
#> 
#> $sources$HP$version
#> [1] "annotations: 01.2024\nclasses: None"
#> 
#> 
#> $sources$HPA
#> $sources$HPA$name
#> [1] "Human Protein Atlas"
#> 
#> $sources$HPA$version
#> [1] "annotations: HPA website: 23-07-17\nclasses: script: 24-01-02"
#> 
#> 
#> $sources$KEGG
#> $sources$KEGG$name
#> [1] "Kyoto Encyclopedia of Genes and Genomes"
#> 
#> $sources$KEGG$version
#> [1] "KEGG FTP Release 2024-01-22"
#> 
#> 
#> $sources$MIRNA
#> $sources$MIRNA$name
#> [1] "miRTarBase"
#> 
#> $sources$MIRNA$version
#> [1] "Release 9.0"
#> 
#> 
#> $sources$REAC
#> $sources$REAC$name
#> [1] "Reactome"
#> 
#> $sources$REAC$version
#> [1] "annotations: BioMart\nclasses: 2024-1-25"
#> 
#> 
#> $sources$TF
#> $sources$TF$name
#> [1] "Transfac"
#> 
#> $sources$TF$version
#> [1] "annotations: TRANSFAC Release 2023.2\nclasses: v2"
#> 
#> 
#> $sources$WP
#> $sources$WP$name
#> [1] "WikiPathways"
#> 
#> $sources$WP$version
#> [1] "20240101"
#> 
#> 
#> 
#> $taxonomy_id
#> [1] "9606"

Custom data sources with upload_GMT_file

In addition to the available GO, KEGG, etc data sources, users can upload their own custom data source using the Gene Matrix Transposed file format (GMT). The file format is described in here. The users can compose the files themselves or use pre-compiled gene sets from available dedicated websites like Molecular Signatures Database (MSigDB), etc. The GMT files for g:Profiler default sources (except for KEGG and Transfac as we are restricted by data source licenses) are downloadabale from the Data sources section in g:Profiler.

upload_GMT_file enables to upload GMT file(s). The input gmtfile is the filename of the GMT file together with the path to the file. The input can also be several GMT files compressed into a ZIP file. The file extension should be .gmt or .zip in case of multiple GMT files. The uploaded filename is used to define the source name in the enrichment results.

For example, using the BioCarta gene sets downloaded from the MSigDB Collections.

download.file(url = "http://software.broadinstitute.org/gsea/resources/msigdb/7.0/c2.cp.biocarta.v7.0.symbols.gmt", destfile = "extdata/biocarta.gmt")

upload_GMT_file(gmtfile = "extdata/biocarta.gmt")

The result is a string that denotes the unique ID of the uploaded data source in the g:Profiler database. In this examaple, the ID is gp_ _TEXF_hZLM_d18.

After the upload, this ID can be used as a value for the parameter organism in the gost function. The input query should consist of identifiers that are available in the GMT file. Note that all the genes in the GMT file define the domain size and therefore it is not sufficient to include only the selection of interesting terms to the file.

custom_gostres <- gost(query = c("MAPK3",   "PIK3C2G", "HRAS", "PIK3R1", "MAP2K1", 
                                 "RAF1", "PLCG1",   "GNAQ", "MAPK1", "PRKCB",   "CRK", "BCAR1", "NFKB1"),
                       organism = "gp__TEXF_hZLM_d18")
#> Detected custom GMT source request
head(custom_gostres$result, 3)
#>     query significant      p_value term_size query_size intersection_size
#> 1 query_1        TRUE 5.324995e-26        19         13                13
#> 2 query_1        TRUE 1.690996e-12        27         13                 9
#> 3 query_1        TRUE 8.094988e-12        19         13                 8
#>   precision    recall                term_id   source
#> 1 1.0000000 0.6842105 BIOCARTA_CXCR4_PATHWAY biocarta
#> 2 0.6923077 0.3333333  BIOCARTA_PYK2_PATHWAY biocarta
#> 3 0.6153846 0.4210526  BIOCARTA_CCR3_PATHWAY biocarta
#>                                                             term_name
#> 1 http://www.gsea-msigdb.org/gsea/msigdb/cards/BIOCARTA_CXCR4_PATHWAY
#> 2  http://www.gsea-msigdb.org/gsea/msigdb/cards/BIOCARTA_PYK2_PATHWAY
#> 3  http://www.gsea-msigdb.org/gsea/msigdb/cards/BIOCARTA_CCR3_PATHWAY
#>   effective_domain_size source_order parents
#> 1                  1476           47    NULL
#> 2                  1476          109    NULL
#> 3                  1476           31    NULL

There is no need to repeatedly upload the same GMT file(s) every time before the enrichment analysis. This can only be uploaded once and then the ID can be used in any further enrichment analyses that are based on that custom source. The same ID can also be used in the web tool as a token under the Custom GMT options. For example, the same query in the web tool is available from https://biit.cs.ut.ee/gplink/l/jh3HdbUWQZ.