Making sense of gene and proteins lists with functional enrichment analysis

10 gprofiler2

gprofiler2 is the R interface to the g:Profiler web-based toolset that you used in day 1 of the workshop.

Like the web interface, gprofiler2 performs ORA with g:GOSt against multiple databases simultaneously.

It supports all the same organisms, namespaces and data sources as the web tool. The list of organisms and corresponding data sources is available here (n = 984).

The full list of supported namespaces is available here.

The gprofiler2 user guide can be found here.

 

10.1 Input data

Since we are doing ORA, we will need a filtered gene list, and a background gene list. We will continue with the RNAseq dataset from Pezzini et al 2017 introduced yesterday.

 

10.2 Activity overview

  1. Load input dataset (a gene matrix with adjusted P values and log2 fold change values)
  2. Filter for differentially expressed genes (DEGs) and create a gene list R object
  3. Extract background genes and create a background gene list R object
  4. Run ORA with gost function
  5. Save the tabular results to a file
  6. Visualise the results with gostplot
  7. Run a gost multi-query for up-regulated and down-regulated genes
  8. Compare gprofiler2 R results to the g:Profiler web results

 

10.3 Analysis workflow

The complete executable analysis is included below. Run the code chunks interactively while working through the chapter, or render the full Bookdown project to execute the analysis and retain the code, parameters, results, figures, and software-version information in the book output.

10.4 0. Working directory

Bookdown evaluates relative file paths from the book project directory. Confirm that the input files used below are available from that directory, or update the file paths to match your project structure:

#> [1] "/home/runner/work/Functional_Enrichment_Web_and_Command_Line/Functional_Enrichment_Web_and_Command_Line"
#> [1] "/home/runner/work/Functional_Enrichment_Web_and_Command_Line/Functional_Enrichment_Web_and_Command_Line"

10.5 1. Load input data

Raw data from Pezzini et al 2017 was subjected to differential gene expression analysis with Degust and the results file saved to Pezzini_DE.txt.

The input data file is within the current working directory so we do not need to specify its directory path.

#> # A tibble: 6 × 10
#>   Gene.ID      Gene.Name  Log2FC     FDR `undif-r1` `undif-r2` `undif-r3` `dif-r1` `dif-r2` `dif-r3`
#>   <chr>        <chr>       <dbl>   <dbl>      <dbl>      <dbl>      <dbl>    <dbl>    <dbl>    <dbl>
#> 1 ENSG0000000… TSPAN6     0.929  9.00e-4        370        343        367      722      774      669
#> 2 ENSG0000000… DPM1       0.0858 4.53e-1       1001        912        998     1004     1136     1107
#> 3 ENSG0000000… SCYL3     -0.0244 8.13e-1        453        421        409      390      463      476
#> 4 ENSG0000000… C1orf112  -1.21   3.78e-3        674        589        603      159      380      353
#> 5 ENSG0000000… CFH        6.01   3.41e-3          0          0          0       70       28       21
#> 6 ENSG0000000… FUCA2      0.730  7.04e-2        499        469        416     1185      690      651

The dataframe shows genes with fold change and FDR values, along with some normalised counts values for the 6 samples (2 groups with 3 replicates each).

Look on the environment pane of RStudio, and you can see a description ‘14420 obs. of 10 variables’ - this shows your dataframe consists of 10 columns and 14,420 genes.

10.6 2. Get ORA gene list

Now we need to filter for differentially expressed genes (DEGs), and we will apply the thresholds adjusted P values/FDR < 0.01, and log2fold change of 2.

We will use ENSEMBL gene IDs (column 1).

#> Number of genes passing FDR and fold change filter: 792
#> Table saved to Pezzini_DEGs.txt

We have 792 genes passing our filters.

10.7 3. Get background gene list

Recall from the webinar and day 1 of the workshop that an experimental background gene list is crucial to avoiding false positives and minimising tissue bias with ORA.

The analysis in Degust has already removed lowly expressed genes, so we can simply extract all genes from this data matrix as our background gene list and save it as our ‘background’ object, as well as save to disk so that we can include it within the supplementary materials of any resultant publications for reproducibility.

#> Number of background genes: 14420
#> Table saved to Pezzini_background_genes.txt

10.8 4. Run ORA with gost function

Before running the below code chunk, review the parameters for the gost ORA function. We can do this easily by bringing up the help page for the function in the Help pane of RStudio. This is the same information that is displayed in the gprofiler2 user guide.

Observe the similarities to the parameters available on the g:Profiler web interface, for example organism, the correction method (g:Profiler’s custom g_scs method), and domain scope (background genes).

Run the below code which explicitly includes all available gost parameters. Including all parameters, even if the defaults suit your needs, makes your parameter choices explicit. Sometimes, default settings can change between versions!

An error-free gost run should produce no console output. As the code is running, there wll be a green bar to the elft of the code chunk.

Our results are saved in the R object ora.

View the top-most significant enrichments with the R head command. Only significant enrichments passing your specified threshold (adjusted P value < 0.05) are included in the results object because we have included significant = TRUE.

Use the black arrow on the right of the table to scroll to other columns.

#> [1] "result" "meta"
#>     query significant      p_value term_size query_size intersection_size precision     recall
#> 1 query_1        TRUE 8.665384e-24      5136        726               430 0.5922865 0.08372274
#> 2 query_1        TRUE 5.229851e-20      3195        726               300 0.4132231 0.09389671
#> 3 query_1        TRUE 1.247572e-17      3706        726               326 0.4490358 0.08796546
#> 4 query_1        TRUE 3.718003e-16      2106        726               215 0.2961433 0.10208927
#> 5 query_1        TRUE 5.767355e-16      4540        726               372 0.5123967 0.08193833
#> 6 query_1        TRUE 1.221690e-15       962        726               126 0.1735537 0.13097713
#>      term_id source                          term_name effective_domain_size source_order
#> 1 GO:0032501  GO:BP   multicellular organismal process                 13009         7116
#> 2 GO:0048731  GO:BP                 system development                 13009        12138
#> 3 GO:0007275  GO:BP multicellular organism development                 13009         2666
#> 4 GO:0009653  GO:BP anatomical structure morphogenesis                 13009         3414
#> 5 GO:0048856  GO:BP   anatomical structure development                 13009        12247
#> 6 GO:0007267  GO:BP                cell-cell signaling                 13009         2659
#>                  parents highlighted
#> 1             GO:0008150        TRUE
#> 2 GO:0007275, GO:0048856       FALSE
#> 3 GO:0032501, GO:0048856       FALSE
#> 4 GO:0032502, GO:0048856       FALSE
#> 5             GO:0032502       FALSE
#> 6 GO:0007154, GO:0023052       FALSE

Let’s give our query a name:

#>               query significant      p_value term_size query_size intersection_size precision
#> 1 DEGs_Padj0.05_FC2        TRUE 8.665384e-24      5136        726               430 0.5922865
#> 2 DEGs_Padj0.05_FC2        TRUE 5.229851e-20      3195        726               300 0.4132231
#> 3 DEGs_Padj0.05_FC2        TRUE 1.247572e-17      3706        726               326 0.4490358
#> 4 DEGs_Padj0.05_FC2        TRUE 3.718003e-16      2106        726               215 0.2961433
#> 5 DEGs_Padj0.05_FC2        TRUE 5.767355e-16      4540        726               372 0.5123967
#> 6 DEGs_Padj0.05_FC2        TRUE 1.221690e-15       962        726               126 0.1735537
#>       recall    term_id source                          term_name effective_domain_size
#> 1 0.08372274 GO:0032501  GO:BP   multicellular organismal process                 13009
#> 2 0.09389671 GO:0048731  GO:BP                 system development                 13009
#> 3 0.08796546 GO:0007275  GO:BP multicellular organism development                 13009
#> 4 0.10208927 GO:0009653  GO:BP anatomical structure morphogenesis                 13009
#> 5 0.08193833 GO:0048856  GO:BP   anatomical structure development                 13009
#> 6 0.13097713 GO:0007267  GO:BP                cell-cell signaling                 13009
#>   source_order                parents highlighted
#> 1         7116             GO:0008150        TRUE
#> 2        12138 GO:0007275, GO:0048856       FALSE
#> 3         2666 GO:0032501, GO:0048856       FALSE
#> 4         3414 GO:0032502, GO:0048856       FALSE
#> 5        12247             GO:0032502       FALSE
#> 6         2659 GO:0007154, GO:0023052       FALSE

We can obtain a list of queried databases:

#> 
#> GO:BP GO:CC GO:MF    HP   HPA  KEGG MIRNA  REAC    TF    WP 
#>   163    22    26     1     3     3     1    16    13     6

Same as the web tool, we have enrichment results for GP (BP, CC, MF), HP (human phenotype), HPA (human protein atlas), KEGG, MiRNA, Reactome, Transcription Factors, and WikiPathways.

10.9 5. Save the results to a file

Let’s save the results file to disk. This is handy when you want to export results elsewhere for further analysis, for example Excel.

First, we will re-order the columns so the output more closely matches the tables that are downloaded from the web version of the tool.

#>   source                          term_name    term_id      p_value term_size query_size
#> 1  GO:BP   multicellular organismal process GO:0032501 8.665384e-24      5136        726
#> 2  GO:BP                 system development GO:0048731 5.229851e-20      3195        726
#> 3  GO:BP multicellular organism development GO:0007275 1.247572e-17      3706        726
#> 4  GO:BP anatomical structure morphogenesis GO:0009653 3.718003e-16      2106        726
#> 5  GO:BP   anatomical structure development GO:0048856 5.767355e-16      4540        726
#> 6  GO:BP                cell-cell signaling GO:0007267 1.221690e-15       962        726
#>   intersection_size effective_domain_size
#> 1               430                 13009
#> 2               300                 13009
#> 3               326                 13009
#> 4               215                 13009
#> 5               372                 13009
#> 6               126                 13009
#> Table saved to gprofiler_ORA_results.csv

gprofiler has a function to print tables that mimic the web tool called publish_gosttable. These are image files, so not for importing to Excel like the CSV we just created.

Let’s extract the results for Reactome and save to a gosttable.

#> The image is saved to gprofiler_Reactome_gosttable.pdf

10.10 6. Visualise the results

10.10.1 Manhattan plots

The gostplot function creates a Manhattan plot similar to the one shown on the web tool. By applying the parameter interactive=TRUE we can hover over the data points to see enriched term details.

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 since p-values smaller than that can all be summarised as ‘highly significant’.

#> Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
#> ℹ Please use the `linewidth` argument instead.
#> ℹ The deprecated feature was likely used in the gprofiler2 package.
#>   Please report the issue at <https://biit.cs.ut.ee/gprofiler/page/contact>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the gprofiler2 package.
#>   Please report the issue at <https://biit.cs.ut.ee/gprofiler/page/contact>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

There are a lot of significant enrichments for GO biological processes. Many of these are probably terms containing a large number of genes, so not particularly informative. Other R tools have default settings limiting the minimum and maximum number of genes in a geneset to be included in the analysis. Since there is no direct parameter to restrict term size to gostplot, we can filter the ORA results before plotting. Let’s apply a maximum gene set size of 500, and a minimum gene set size of 10, which are the default setting used by clusterProfiler.

This has cleaned up ‘Biological Process’ a little bit, enabling signals of more specific terms to be highlighted.

gprofiler2 includes a function for creating a publication-ready image that can optionally highlight specific terms. We need to first produce a plot with interactice = FALSE, save it to an object, and then provide that plot object to the publish_gostplot function.

The publish_gostplot parameter highlight_terms enables you to highlight specific terms on the plot, with a table showing enrichment details below for those highlighted terms.

Let’s highlight some selected terms manually. You need to provide the term ID not term name.

#> The image is saved to gprofiler_collagen_gostplot.pdf

Like g:Profiler web, the coloured boxes on the table are by adjusted P value, with darker colours indicating more significant results. Colours range from yellow through green to dark blue.

You can use R grepl function to search for terms with names matching some keyword. Let’s highlight all terms related to receptors. The code chunk applies an increased figure height, to ensure we can see the whole plot within the rendered chapter.

#> The image is saved to gprofiler_receptors_gostplot.pdf

10.10.2 Dotplots

One of the advantages of working in R is flexibility with visualisations. While the interactive Manhattan plots and publish_gostplot options are nice, it can also be useful to visualise P values against all term descriptions.

One way to do this is with a dotplot. We can loop through all databases and use the R package ggplot2 to make a dotplot for each database with significantly enriched terms for our gene list.

#> No significant enrichments for database: HP

#> No significant enrichments for database: MIRNA
#> No significant enrichments for database: TF

For plots with a lot of enriched terms, such as GO Biological Process, the display within the rendered chapter is less than ideal. Saving the plot to an image file enables better resolution:

#> png 
#>   2
#> Table saved to gprofiler_GO_BP_dotplot.pdf

Open this plot by clicking it from the ‘Files’ pane of RStudio. Notice how the term names are now readable :-)

10.11 7. Run a gost multi-query

By providing more than one gene list and setting multi_query = TRUE, results from all of the gene lists are grouped by term IDs for easier comparison. This can be handy when you have multiple comparisons within an experiment, or when you want to investigate enrichments within the up and down regulated genes separately.

First, we need to extracts separate gene lists for up-regulated and down-regulated genes.

#> Number of upregulated DEGs: 577
#> Up-regulated DEGs saved to up_DEGs.txt
#> Number of downregulated DEGs: 215
#> Down-regulated DEGs saved to down_DEGs.txt

Now run gost as multi-query. This may take a few moments.

The changes required for multi-query are providing a list of gene list objects to the query parameter instead of a single gene list object, and setting multi_query = TRUE, which is FALSE by default. By including all genes as well, we can efficiently compare up vs down vs no separation.

Now create a multi-query interactive Manhattan plot with gostplot:

Unfortunately, the rendered chapter squashes the top plot over the bottom one, and adjusting figure height or plot layout options doesn’t seem to help. Plotting as non-interactive or plotting from the console to the plots pane both produce a correct looking plot.

#> The image is saved to gprofiler_ORA_multiquery.pdf

To access the tabular results separately, they need to be split, as a number of the columns are comma-delimited lists with one value for each of the 3 queries.

#>      term_id                                 p_values       significant term_size   query_sizes
#> 1 GO:0005576 6.964351e-25, 1.000000e+00, 1.978590e-28 TRUE, FALSE, TRUE      2726 559, 205, 764
#> 2 GO:0071944 5.417977e-18, 3.190965e-03, 7.003641e-25  TRUE, TRUE, TRUE      4001 559, 205, 764
#> 3 GO:0032501 5.800468e-13, 1.208643e-07, 8.665384e-24  TRUE, TRUE, TRUE      5136 529, 197, 726
#> 4 GO:0005615 3.451654e-21, 1.000000e+00, 1.390304e-22 TRUE, FALSE, TRUE      2225 559, 205, 764
#> 5 GO:0031012 3.739372e-22, 1.000000e+00, 1.693574e-22 TRUE, FALSE, TRUE       244 559, 205, 764
#> 6 GO:0030312 3.739372e-22, 1.000000e+00, 1.693574e-22 TRUE, FALSE, TRUE       244 559, 205, 764
#>   intersection_sizes source                        term_name effective_domain_size source_order
#> 1       225, 66, 291  GO:CC             extracellular region                 13682          194
#> 2       271, 96, 367  GO:CC                   cell periphery                 13682         2724
#> 3      304, 126, 430  GO:BP multicellular organismal process                 13009         7116
#> 4       190, 51, 241  GO:CC              extracellular space                 13682          227
#> 5          56, 9, 65  GO:CC             extracellular matrix                 13682         1000
#> 6          56, 9, 65  GO:CC external encapsulating structure                 13682          902
#>                  parents
#> 1             GO:0110165
#> 2             GO:0110165
#> 3             GO:0008150
#> 4 GO:0005576, GO:0110165
#> 5             GO:0030312
#> 6 GO:0071944, GO:0110165

If you run the command head(ora_multi$result) directly in the console (not from the chapter) you can see the list values.

One might want to explore the comparisons in more detail for example viewing the separate P values in Excel, so having these results as a file would be handy. The columns that are lists need to be converted to character format before printing to TSV.

#>               term_id              p_values           significant             term_size 
#>           "character"                "list"                "list"             "integer" 
#>           query_sizes    intersection_sizes                source             term_name 
#>                "list"                "list"           "character"           "character" 
#> effective_domain_size          source_order               parents 
#>             "integer"             "integer"                "list"

Convert the columns that are lists to characters

#>               term_id              p_values           significant             term_size 
#>           "character"           "character"           "character"             "integer" 
#>           query_sizes    intersection_sizes                source             term_name 
#>           "character"           "character"           "character"           "character" 
#> effective_domain_size          source_order               parents 
#>             "integer"             "integer"           "character"

View the new table format:

#>      term_id                                                      p_values     significant
#> 1 GO:0005576                   6.96435117941893e-25,1,1.97858995909489e-28 TRUE,FALSE,TRUE
#> 2 GO:0071944 5.41797719942088e-18,0.00319096483707521,7.00364071713173e-25  TRUE,TRUE,TRUE
#> 3 GO:0032501 5.8004684270304e-13,1.20864284962376e-07,8.66538362719036e-24  TRUE,TRUE,TRUE
#> 4 GO:0005615                   3.45165351553507e-21,1,1.39030399342413e-22 TRUE,FALSE,TRUE
#> 5 GO:0031012                   3.73937179374249e-22,1,1.69357410889906e-22 TRUE,FALSE,TRUE
#> 6 GO:0030312                   3.73937179374249e-22,1,1.69357410889906e-22 TRUE,FALSE,TRUE
#>   term_size query_sizes intersection_sizes source                        term_name
#> 1      2726 559,205,764         225,66,291  GO:CC             extracellular region
#> 2      4001 559,205,764         271,96,367  GO:CC                   cell periphery
#> 3      5136 529,197,726        304,126,430  GO:BP multicellular organismal process
#> 4      2225 559,205,764         190,51,241  GO:CC              extracellular space
#> 5       244 559,205,764            56,9,65  GO:CC             extracellular matrix
#> 6       244 559,205,764            56,9,65  GO:CC external encapsulating structure
#>   effective_domain_size source_order               parents
#> 1                 13682          194            GO:0110165
#> 2                 13682         2724            GO:0110165
#> 3                 13009         7116            GO:0008150
#> 4                 13682          227 GO:0005576,GO:0110165
#> 5                 13682         1000            GO:0030312
#> 6                 13682          902 GO:0071944,GO:0110165
#> ORA multi-query results written to gprofiler_ORA_ora_multiquery.tsv

10.12 8. Compare gprofiler2 R results to g:Profiler web results

In day 1 of the workshop, you ran ORA with g:Profiler web tool and saved the results to a CSV. Let’s compare the results to those we have generated in R. Do we expect the results to be identical or differ slightly?

The input file here is one that we have created, but should match yours as long as you used the same P filters and gprofiler parameters.

#>   source                                 term_name    term_id highlighted adjusted_p_value
#> 1  GO:MF               signaling receptor activity GO:0038023        true     4.798199e-16
#> 2  GO:MF             molecular transducer activity GO:0060089       false     4.798199e-16
#> 3  GO:MF transmembrane signaling receptor activity GO:0004888       false     1.143972e-14
#> 4  GO:MF                 glycosaminoglycan binding GO:0005539        true     1.559822e-11
#> 5  GO:MF     signaling receptor activator activity GO:0030546        true     4.269590e-11
#> 6  GO:MF     signaling receptor regulator activity GO:0030545       false     7.017863e-11
#>   negative_log10_of_adjusted_p_value term_size query_size intersection_size effective_domain_size
#> 1                           15.31892       626        748                96                 13659
#> 2                           15.31892       626        748                96                 13659
#> 3                           13.94158       474        748                79                 13659
#> 4                           10.80693       154        748                39                 13659
#> 5                           10.36961       242        748                49                 13659
#> 6                           10.15380       263        748                51                 13659
#>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     intersections
#> 1 ENSG00000006016,ENSG00000006071,ENSG00000011028,ENSG00000013588,ENSG00000021645,ENSG00000049249,ENSG00000049323,ENSG00000064300,ENSG00000064989,ENSG00000066468,ENSG00000067182,ENSG00000069431,ENSG00000075275,ENSG00000077092,ENSG00000077238,ENSG00000083454,ENSG00000091592,ENSG00000101000,ENSG00000101292,ENSG00000104213,ENSG00000105464,ENSG00000106991,ENSG00000110195,ENSG00000111199,ENSG00000112038,ENSG00000112164,ENSG00000113763,ENSG00000115353,ENSG00000116678,ENSG00000117525,ENSG00000119508,ENSG00000121966,ENSG00000122679,ENSG00000123146,ENSG00000125384,ENSG00000127418,ENSG00000128052,ENSG00000132170,ENSG00000132256,ENSG00000135902,ENSG00000136160,ENSG00000144476,ENSG00000145451,ENSG00000145623,ENSG00000148053,ENSG00000150594,ENSG00000152894,ENSG00000153029,ENSG00000153294,ENSG00000153707,ENSG00000154133,ENSG00000154928,ENSG00000156687,ENSG00000157404,ENSG00000158270,ENSG00000160013,ENSG00000162894,ENSG00000163394,ENSG00000163513,ENSG00000164082,ENSG00000164342,ENSG00000164509,ENSG00000166206,ENSG00000166825,ENSG00000167664,ENSG00000168412,ENSG00000168539,ENSG00000169860,ENSG00000170775,ENSG00000171631,ENSG00000172572,ENSG00000173567,ENSG00000173662,ENSG00000176884,ENSG00000179546,ENSG00000180720,ENSG00000180739,ENSG00000180772,ENSG00000181072,ENSG00000181656,ENSG00000182771,ENSG00000183098,ENSG00000185291,ENSG00000185551,ENSG00000185924,ENSG00000186827,ENSG00000187323,ENSG00000196187,ENSG00000196277,ENSG00000196664,ENSG00000198121,ENSG00000198963,ENSG00000213694,ENSG00000221866,ENSG00000237988,ENSG00000243364
#> 2 ENSG00000006016,ENSG00000006071,ENSG00000011028,ENSG00000013588,ENSG00000021645,ENSG00000049249,ENSG00000049323,ENSG00000064300,ENSG00000064989,ENSG00000066468,ENSG00000067182,ENSG00000069431,ENSG00000075275,ENSG00000077092,ENSG00000077238,ENSG00000083454,ENSG00000091592,ENSG00000101000,ENSG00000101292,ENSG00000104213,ENSG00000105464,ENSG00000106991,ENSG00000110195,ENSG00000111199,ENSG00000112038,ENSG00000112164,ENSG00000113763,ENSG00000115353,ENSG00000116678,ENSG00000117525,ENSG00000119508,ENSG00000121966,ENSG00000122679,ENSG00000123146,ENSG00000125384,ENSG00000127418,ENSG00000128052,ENSG00000132170,ENSG00000132256,ENSG00000135902,ENSG00000136160,ENSG00000144476,ENSG00000145451,ENSG00000145623,ENSG00000148053,ENSG00000150594,ENSG00000152894,ENSG00000153029,ENSG00000153294,ENSG00000153707,ENSG00000154133,ENSG00000154928,ENSG00000156687,ENSG00000157404,ENSG00000158270,ENSG00000160013,ENSG00000162894,ENSG00000163394,ENSG00000163513,ENSG00000164082,ENSG00000164342,ENSG00000164509,ENSG00000166206,ENSG00000166825,ENSG00000167664,ENSG00000168412,ENSG00000168539,ENSG00000169860,ENSG00000170775,ENSG00000171631,ENSG00000172572,ENSG00000173567,ENSG00000173662,ENSG00000176884,ENSG00000179546,ENSG00000180720,ENSG00000180739,ENSG00000180772,ENSG00000181072,ENSG00000181656,ENSG00000182771,ENSG00000183098,ENSG00000185291,ENSG00000185551,ENSG00000185924,ENSG00000186827,ENSG00000187323,ENSG00000196187,ENSG00000196277,ENSG00000196664,ENSG00000198121,ENSG00000198963,ENSG00000213694,ENSG00000221866,ENSG00000237988,ENSG00000243364
#> 3                                                                                                                                                                                                                                                                                 ENSG00000006016,ENSG00000006071,ENSG00000013588,ENSG00000021645,ENSG00000049323,ENSG00000064300,ENSG00000064989,ENSG00000066468,ENSG00000067182,ENSG00000069431,ENSG00000075275,ENSG00000077238,ENSG00000083454,ENSG00000101292,ENSG00000104213,ENSG00000105464,ENSG00000106991,ENSG00000111199,ENSG00000112038,ENSG00000112164,ENSG00000113763,ENSG00000115353,ENSG00000116678,ENSG00000117525,ENSG00000121966,ENSG00000122679,ENSG00000123146,ENSG00000125384,ENSG00000127418,ENSG00000128052,ENSG00000132170,ENSG00000135902,ENSG00000136160,ENSG00000144476,ENSG00000145451,ENSG00000145623,ENSG00000148053,ENSG00000150594,ENSG00000152894,ENSG00000153029,ENSG00000153294,ENSG00000153707,ENSG00000154928,ENSG00000156687,ENSG00000157404,ENSG00000160013,ENSG00000162894,ENSG00000163394,ENSG00000163513,ENSG00000164082,ENSG00000164342,ENSG00000164509,ENSG00000166206,ENSG00000168412,ENSG00000168539,ENSG00000169860,ENSG00000170775,ENSG00000171631,ENSG00000173567,ENSG00000173662,ENSG00000176884,ENSG00000179546,ENSG00000180720,ENSG00000180739,ENSG00000180772,ENSG00000181072,ENSG00000181656,ENSG00000182771,ENSG00000185291,ENSG00000186827,ENSG00000187323,ENSG00000196187,ENSG00000196277,ENSG00000198121,ENSG00000198963,ENSG00000213694,ENSG00000221866,ENSG00000237988,ENSG00000243364
#> 4                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 ENSG00000000971,ENSG00000009694,ENSG00000011465,ENSG00000012223,ENSG00000038427,ENSG00000066468,ENSG00000080573,ENSG00000091986,ENSG00000101282,ENSG00000104332,ENSG00000106991,ENSG00000108688,ENSG00000108700,ENSG00000110492,ENSG00000115414,ENSG00000123610,ENSG00000127418,ENSG00000130203,ENSG00000130635,ENSG00000135919,ENSG00000140285,ENSG00000154175,ENSG00000163430,ENSG00000163513,ENSG00000165973,ENSG00000166106,ENSG00000169245,ENSG00000169248,ENSG00000169429,ENSG00000172638,ENSG00000173376,ENSG00000180875,ENSG00000182492,ENSG00000184613,ENSG00000185924,ENSG00000188783,ENSG00000196562,ENSG00000204248,ENSG00000204381
#> 5                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 ENSG00000001617,ENSG00000006016,ENSG00000019991,ENSG00000075213,ENSG00000078401,ENSG00000078579,ENSG00000081041,ENSG00000092421,ENSG00000102265,ENSG00000104826,ENSG00000105329,ENSG00000106991,ENSG00000108688,ENSG00000108691,ENSG00000108700,ENSG00000110492,ENSG00000113739,ENSG00000114251,ENSG00000119699,ENSG00000120659,ENSG00000125538,ENSG00000125848,ENSG00000137077,ENSG00000139269,ENSG00000140285,ENSG00000143816,ENSG00000143869,ENSG00000147869,ENSG00000150281,ENSG00000153993,ENSG00000162344,ENSG00000163273,ENSG00000163377,ENSG00000164326,ENSG00000165197,ENSG00000166923,ENSG00000169245,ENSG00000169248,ENSG00000169429,ENSG00000172247,ENSG00000175505,ENSG00000180875,ENSG00000189292,ENSG00000196154,ENSG00000198719,ENSG00000219438,ENSG00000232810,ENSG00000266524,ENSG00000271503
#> 6                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 ENSG00000001617,ENSG00000006016,ENSG00000019991,ENSG00000075213,ENSG00000078401,ENSG00000078579,ENSG00000081041,ENSG00000092421,ENSG00000102265,ENSG00000104826,ENSG00000105329,ENSG00000106991,ENSG00000108688,ENSG00000108691,ENSG00000108700,ENSG00000110492,ENSG00000113739,ENSG00000114251,ENSG00000119699,ENSG00000120659,ENSG00000125538,ENSG00000125848,ENSG00000137077,ENSG00000139269,ENSG00000140285,ENSG00000143816,ENSG00000143869,ENSG00000147255,ENSG00000147869,ENSG00000150281,ENSG00000153993,ENSG00000162344,ENSG00000163273,ENSG00000163377,ENSG00000164326,ENSG00000165197,ENSG00000166923,ENSG00000169245,ENSG00000169248,ENSG00000169429,ENSG00000172247,ENSG00000175505,ENSG00000180772,ENSG00000180875,ENSG00000189292,ENSG00000196154,ENSG00000198719,ENSG00000219438,ENSG00000232810,ENSG00000266524,ENSG00000271503

Check the numbers: are there any terms significant from one tool but not the other?

#> [1] "Number of significant terms from web: 273"
#> [1] "Number of significant terms from R: 254"
#> [1] "Number of terms unique to web: 45"
#> [1] "Number of terms unique to gprofiler2 (R): 26"

That’s a good start! Do the P values differ? Let’s look closely at the GO ‘Molecular Function’ P values via a barplot.

Format the P values for plotting:

Create barplot to compare P values web vs R:

Great! We know we applied the same parameters, and used the same input gene lists. The identical results must mean that gprofiler2 is using the same database version as g:Profiler web.

10.13 9. Save versions and session details

10.13.1 gProfiler database version

Let’s check: yesterday when you ran ORA on the web, hopefully you saved your ‘query parameters’ as well as your results.

From my run, I can see version as ‘e111_eg58_p18_f463989d’.

Let’s report the g:Profiler database version used in our analysis:

#> [1] "g:Profiler database version: e114_eg62_p19_27110d83"
#> [1] "gprofiler2 package version: 0.2.4"

10.13.2 R version and R package versions

We can also capture the version of R and other session details including all loaded packages and versions with the sessionInfo() function:

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8    
#>  [5] LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8       LC_NAME=C             
#>  [9] LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices datasets  utils     methods   base     
#> 
#> other attached packages:
#>  [1] tidyr_1.3.2      ggplot2_4.0.3    gprofiler2_0.2.4 readr_2.2.0      lubridate_1.9.5 
#>  [6] dplyr_1.2.1      downlit_0.4.5    kableExtra_1.4.1 knitr_1.51       fs_2.1.0        
#> [11] yaml_2.3.12     
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        xfun_0.60           bslib_0.11.0        htmlwidgets_1.6.4  
#>  [5] tzdb_0.5.0          crosstalk_1.2.2     bitops_1.0-9        vctrs_0.7.3        
#>  [9] tools_4.6.1         generics_0.1.4      curl_7.1.0          parallel_4.6.1     
#> [13] tibble_3.3.1        pkgconfig_2.0.3     data.table_1.18.4   RColorBrewer_1.1-3 
#> [17] S7_0.2.2            lifecycle_1.0.5     compiler_4.6.1      farver_2.1.2       
#> [21] stringr_1.6.0       textshaping_1.0.5   htmltools_0.5.9     sass_0.4.10        
#> [25] lazyeval_0.2.3      RCurl_1.98-1.19     plotly_4.12.1       pillar_1.11.1      
#> [29] crayon_1.5.3        jquerylib_0.1.4     cachem_1.1.0        tidyselect_1.2.1   
#> [33] digest_0.6.39       stringi_1.8.7       purrr_1.2.2         bookdown_0.47      
#> [37] labeling_0.4.3      fastmap_1.2.0       grid_4.6.1          cli_3.6.6          
#> [41] magrittr_2.0.5      utf8_1.2.6          withr_3.0.3         scales_1.4.0       
#> [45] bit64_4.8.2         timechange_0.4.0    rmarkdown_2.31      httr_1.4.8         
#> [49] bit_4.6.0           otel_0.2.0          gridExtra_2.3.1     png_0.1-9          
#> [53] hms_1.1.4           memoise_2.0.1       evaluate_1.0.5      viridisLite_0.4.3  
#> [57] rlang_1.3.0         glue_1.8.1          BiocManager_1.30.27 xml2_1.6.0         
#> [61] renv_1.1.5          svglite_2.2.2       rstudioapi_0.19.0   vroom_1.7.1        
#> [65] jsonlite_2.0.0      R6_2.6.1            systemfonts_1.3.2

10.13.3 RStudio version

Typically, we would simply run RStudio.Version() to print the version details. However, when we knit this document to HTML, the RStudio.Version() function is not available and will cause an error. So to make sure our version details are saved to our static record of the work, we will save to a file, then print the file contents back into the chapter.

10.14 Build the Bookdown output

Save all changes before rendering. Build the complete book from the RStudio Build pane using Build Book, or run the project-specific Bookdown render command used by this workshop.

The chapter will render successfully only when all required packages and input files are available and the executable code chunks complete without errors. Output tables and figures written by the analysis will be saved relative to the Bookdown project directory unless a different path is specified.

❗ If the build fails, use the first error message in the render log to identify the missing file, package, object, or failed code chunk.

10.15 End of activity summary

  • We have extracted a gene list and background gene list from a DE dataset and run ORA with gprofiler2 gost function
  • We have plotted the data with gostplot Manhattan plots and ggplot2 dotplots
  • We have run a gost multi-query separating up and down regulated genes
  • We have verified that gprofiler2 results match the results from g:Profiler web
  • We have incorporated the complete analysis into the rendered Bookdown output as a static record of all executed code, analysis parameters, results, and version information for R, RStudio, R packages, and the g:Profiler database