Module 8.2 Lab: Defined gene list (g:Profiler)
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 Unported License. This means that you are able to copy, share and modify the work, as long as the result is distributed under the same license.
Authors: Veronique Voisin, Ruth Isserlin
Presenter: Veronique Voisin
Introduction
During this practical lab, we will perform pathway enrichment analysis using a defined gene list.
Goal of the exercise 1
For this exercise, our goal is to run pathway enrichment analysis using g:Profiler and explore the results. g:Profiler performs a gene-set enrichment analysis using a hypergeometric test (Fisher’s exact test). The Gene Ontology Biological Process, Reactome and WikiPathways sources are going to be used as pathway databases. We wil also upload our custom geneset file.
g:Profiler can be used using the website at https://biit.cs.ut.ee/gprofiler/gost. However, for this practical lab, we will run it from the g:Profiler R package.
We will run the query and explore the table of results and visualize the results as bar and dot plots.
One of the greatest features of g:Profiler is that it is updated on a regular basis and most of the previous versions are available online ont the gprofiler archive.
The gprofielr2 -g:Profiler R implementation is a wrapper for the web version. You require an internet connection to get enrichment results.
Data
g:Profiler requires a list of genes. For this analysis, we used genes we found to be differentially expressed in palbociclib-resistant samples compared with parental samples from the breast cancer cell line MCF7 (see chapter 1.
Exercise 1 - prepare the gene lists
This exercise requires the file “edgeR_resistant_vs_parental_all_genes.csv” located in the data folder.
Step 1 - Install and Load libraries
# List of CRAN packages
cran_packages <- c(
"tidyverse",
"knitr",
"kableExtra",
"glue",
"patchwork",
"ggdendro",
"pheatmap",
"gprofiler2"
)
# Install CRAN packages if not already installed
for (pkg in cran_packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
}
# Load the libraries
library(tidyverse)
library(knitr)
library(kableExtra)
library(glue)
library(patchwork)
library(ggdendro)
library(pheatmap)
library(gprofiler2)
# create the results folder (if it does not exist yet) to save the outputs
# dir.create("results", showWarnings = FALSE)Step 3 - Run g:Profiler
The next step is to run pathway enrichment analysis using gost g:Profiler using the g:Profiler2 R package. For detailed descriptions of all the parameters that can be specified for the gost g:profiler function see the R package information - here and here.
For this query we are specifying -
- query - the set of genes of interest, as loaded in from the gene list file (
query_set). - significant - set to TRUE so that g:Profiler returns only the pathways that it deems significant (adjusted p-value below its default threshold of 0.05). We will filter these results further in a next step.
- ordered_query - set to FALSE (because we are not taking into account the order of the list)
- exclude_iea - set to TRUE. We are removing the electronic inferred annotations from the GO database
- correction_method - set to fdr. by default g:Profiler uses g:Scs
- organism - set to “hsapiens” for homo sapiens. Organism names are constructed by concatenating the first letter of the name and the family name (according to gprofiler2 documentation)
- source - the geneset source databases to use for the analysis. We recommend using GO biological process (GO:BP), WikiPathways (WP) and Reactome (Reac) but there are additional sources you can add (GO molecular function or cellular component(GO:MF, GO:CC), KEGG, transcription factors (TF), microRNA targets (MIRNA), corum complexes (CORUM), Human protein atlas (HPA),Human phenotype ontology (HP) )
- evcodes - set to TRUE so that g:Profiler also returns the genes of our list found in each pathway (column
intersection, genes separated by a comma). It makes the query a little slower.
observe the results
#get the gprofiler results table
enrichment_results <- gprofiler_results_UPgenes$result
#display the names of all columns
colnames(enrichment_results)## [1] "query" "significant" "p_value"
## [4] "term_size" "query_size" "intersection_size"
## [7] "precision" "recall" "term_id"
## [10] "source" "term_name" "effective_domain_size"
## [13] "source_order" "parents" "evidence_codes"
## [16] "intersection"## query significant p_value term_size query_size intersection_size
## 1 query_1 TRUE 0.005977645 635 74 13
## 2 query_1 TRUE 0.005977645 509 74 11
## 3 query_1 TRUE 0.005977645 61 74 5
## 4 query_1 TRUE 0.005977645 1013 74 16
## 5 query_1 TRUE 0.005977645 188 74 7
## 6 query_1 TRUE 0.005977645 73 74 5
## precision recall term_id source
## 1 0.17567568 0.02047244 GO:0035295 GO:BP
## 2 0.14864865 0.02161100 GO:0001944 GO:BP
## 3 0.06756757 0.08196721 GO:0050818 GO:BP
## 4 0.21621622 0.01579467 GO:2000026 GO:BP
## 5 0.09459459 0.03723404 GO:0050817 GO:BP
## 6 0.06756757 0.06849315 GO:0032963 GO:BP
## term_name effective_domain_size
## 1 tube development 16180
## 2 vasculature development 16180
## 3 regulation of coagulation 16180
## 4 regulation of multicellular organismal development 16180
## 5 coagulation 16180
## 6 collagen metabolic process 16180
## source_order parents
## 1 8374 GO:0007275, GO:0048856
## 2 542 GO:0048731, GO:0072359
## 3 12418 GO:0050817, GO:0051239
## 4 23538 GO:0007275, GO:0050793, GO:0051239
## 5 12417 GO:0032501
## 6 7434 GO:0008152
## evidence_codes
## 1 ISS,IMP IGI IEP IBA,IDA,IDA,IEP,IMP IBA TAS,ISS,ISS,IMP,ISS,TAS,IDA,IMP
## 2 ISS,IMP IGI IEP IBA,IDA,IDA,IEP,IMP IBA TAS,IMP ISS,ISS,ISS,TAS,IDA
## 3 IMP,IDA IMP IBA TAS NAS IC,IC,IDA IBA,IMP IBA
## 4 ISS,IDA,IDA IMP IGI IBA,IDA,IDA,IMP,IDA,IMP,IMP,IDA,ISS,ISS IBA,IBA,TAS,IDA,IDA IMP IBA
## 5 IMP,IBA,IDA IMP IBA TAS NAS IC,IDA IC,IDA IBA,IDA,IMP IBA
## 6 IDA IBA,ISS IBA,IMP,IDA,IDA
## intersection
## 1 CAV1,SERPINE1,F3,ADGRB1,TGFBI,ITGA5,CCBE1,LOX,WNK4,ERRFI1,ANGPTL4,RAMP1,CLMP
## 2 CAV1,SERPINE1,F3,ADGRB1,TGFBI,ITGA5,CCBE1,LOX,ERRFI1,ANGPTL4,RAMP1
## 3 CAV1,SERPINE1,F3,F2R,PROCR
## 4 CAV1,AXL,SERPINE1,F3,ADGRB1,DPYSL5,ADAMTS7,BIN1,ITGA5,FBN1,L1CAM,ERRFI1,CSF3,ANGPTL4,SEMA7A,TESC
## 5 CAV1,AXL,SERPINE1,F3,F2R,FBLN1,PROCR
## 6 P3H2,P3H3,COL5A1,F2R,MRC2YOUR TURN: run g:Profiler with the down genes (use same parameters as above)
##Tip: remove eval=FALSE chunk option to run the code
query_set = XXXX
gprofiler_results_DOWNgenes <- gost(query = XXXX,
significant = XXXX,
ordered_query = XXXX,
exclude_iea = XXXX,
correction_method = XXXX,
organism = XXXX,
source = c("XXXX","XXXX","XXXX"),
evcodes = XXXX)
##observe the results (same as for the UP genes)Clear visualization of the results
The enrichment_results dataframe contains different columns. Let’s rearrange the table to display the most important information. We will arrange the table so it contains: the database origin, the name of each of the pathway, term size , query size, intersection size and adjusted pvalue. -
- Term size: number of genes in the pathway as it is in the database
- Query size: number of genes in our gene list
- Intersection size: number of genes overlapping between our gene list and the tested pathway
- pvalue: adjusted pvalue (corrected for multiple hypothesis testing using the Benjamini-Hochberg method)
The results are ordered by significance, from the most significant pathway (lowest p-value) to the least significant.
enrichment_results <- gprofiler_results_UPgenes$result
enrichment_results$p_value <- formatC(enrichment_results$p_value, format = "e", digits = 2)
enrichment_results = enrichment_results[ , c("source", "term_name", "term_size", "query_size", "intersection_size", "p_value")]
kable(head(enrichment_results, n=20), caption = "Enrichment Result") %>%
kable_styling(bootstrap_options = c("striped", "condensed"), font_size = 10)| source | term_name | term_size | query_size | intersection_size | p_value |
|---|---|---|---|---|---|
| GO:BP | tube development | 635 | 74 | 13 | 5.98e-03 |
| GO:BP | vasculature development | 509 | 74 | 11 | 5.98e-03 |
| GO:BP | regulation of coagulation | 61 | 74 | 5 | 5.98e-03 |
| GO:BP | regulation of multicellular organismal development | 1013 | 74 | 16 | 5.98e-03 |
| GO:BP | coagulation | 188 | 74 | 7 | 5.98e-03 |
| GO:BP | collagen metabolic process | 73 | 74 | 5 | 5.98e-03 |
| GO:BP | anatomical structure morphogenesis | 1847 | 74 | 22 | 5.98e-03 |
| GO:BP | cell adhesion mediated by integrin | 80 | 74 | 5 | 7.04e-03 |
| GO:BP | blood vessel morphogenesis | 457 | 74 | 10 | 8.11e-03 |
| GO:BP | tube morphogenesis | 562 | 74 | 11 | 8.11e-03 |
| GO:BP | regulation of multicellular organismal process | 2293 | 74 | 24 | 8.15e-03 |
| GO:BP | blood vessel development | 487 | 74 | 10 | 1.04e-02 |
| GO:BP | positive regulation of blood coagulation | 21 | 74 | 3 | 1.44e-02 |
| GO:BP | positive regulation of hemostasis | 21 | 74 | 3 | 1.44e-02 |
| GO:BP | circulatory system development | 745 | 74 | 12 | 1.45e-02 |
| GO:BP | regulation of blood coagulation | 57 | 74 | 4 | 1.45e-02 |
| GO:BP | regulation of hemostasis | 58 | 74 | 4 | 1.46e-02 |
| GO:BP | wound healing | 340 | 74 | 8 | 1.46e-02 |
| GO:BP | positive regulation of coagulation | 23 | 74 | 3 | 1.46e-02 |
| GO:BP | blood coagulation | 182 | 74 | 6 | 1.61e-02 |
YOUR TURN: visualize the results with the down genes (use same parameters as above)
##Tip: remove eval=FALSE chunk option to run the code
enrichment_results <- XXXX
enrichment_results$p_value <- formatC(XXXX)
enrichment_results = enrichment_results[ , c("XXXX", "XXXX", "XXXX", "XXXX", "XXXX", "XXXX")]
kable(head(enrichment_results, n=20), caption = "Enrichment Result") %>%
kable_styling(bootstrap_options = c("striped", "condensed"), font_size = 10)Step 4 - Explore top table of results
For this analysis, we chose to use 3 databases: GO:BP, Reactome and Wikipathways. If we look at the top table only, we are only going to see the results of GO:BP but it is interesting to look at the top results for each of the 3 databunique(enrichment_results$source)
enrichment_results <- gprofiler_results_UPgenes$result
##Get the name of the databases
unique(enrichment_results$source) ## [1] "GO:BP" "REAC" "WP"##Filter the dataframe to retrieve only the Reactome results
enrichment_results_Reac <- enrichment_results %>%
filter(source == "REAC")
##Now display the top 20 pathways
kable(head(enrichment_results_Reac, n=20), caption = "Enrichment Result") %>%
kable_styling(bootstrap_options = c("striped", "condensed"), font_size = 10)| query | significant | p_value | term_size | query_size | intersection_size | precision | recall | term_id | source | term_name | effective_domain_size | source_order | parents | evidence_codes | intersection |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| query_1 | TRUE | 0.0000077 | 298 | 57 | 12 | 0.2105263 | 0.0402685 | REAC:R-HSA-1474244 | REAC | Extracellular matrix organization | 11056 | 930 | REAC:0000000 | REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC | COL6A1,SERPINE1,P3H2,TIMP2,P3H3,COL6A2,COL5A1,ITGA5,FBLN1,FBN1,LOX,BGN |
| query_1 | TRUE | 0.0008208 | 89 | 57 | 6 | 0.1052632 | 0.0674157 | REAC:R-HSA-1474290 | REAC | Collagen formation | 11056 | 403 | REAC:R-HSA-1474244 | REAC,REAC,REAC,REAC,REAC,REAC | COL6A1,P3H2,P3H3,COL6A2,COL5A1,LOX |
| query_1 | TRUE | 0.0021054 | 67 | 57 | 5 | 0.0877193 | 0.0746269 | REAC:R-HSA-1650814 | REAC | Collagen biosynthesis and modifying enzymes | 11056 | 400 | REAC:R-HSA-1474290 | REAC,REAC,REAC,REAC,REAC | COL6A1,P3H2,P3H3,COL6A2,COL5A1 |
| query_1 | TRUE | 0.0029160 | 76 | 57 | 5 | 0.0877193 | 0.0657895 | REAC:R-HSA-3000178 | REAC | ECM proteoglycans | 11056 | 851 | REAC:R-HSA-1474244 | REAC,REAC,REAC,REAC,REAC | COL6A1,SERPINE1,COL6A2,COL5A1,BGN |
| query_1 | TRUE | 0.0033756 | 44 | 57 | 4 | 0.0701754 | 0.0909091 | REAC:R-HSA-1566948 | REAC | Elastic fibre formation | 11056 | 875 | REAC:R-HSA-1474244 | REAC,REAC,REAC,REAC | ITGA5,FBLN1,FBN1,LOX |
| query_1 | TRUE | 0.0033756 | 84 | 57 | 5 | 0.0877193 | 0.0595238 | REAC:R-HSA-216083 | REAC | Integrin cell surface interactions | 11056 | 1267 | REAC:R-HSA-1474244 | REAC,REAC,REAC,REAC,REAC | COL6A1,COL6A2,COL5A1,ITGA5,FBN1 |
| query_1 | TRUE | 0.0097759 | 60 | 57 | 4 | 0.0701754 | 0.0666667 | REAC:R-HSA-2022090 | REAC | Assembly of collagen fibrils and other multimeric structures | 11056 | 186 | REAC:R-HSA-1474290 | REAC,REAC,REAC,REAC | COL6A1,COL6A2,COL5A1,LOX |
| query_1 | TRUE | 0.0256174 | 140 | 57 | 5 | 0.0877193 | 0.0357143 | REAC:R-HSA-1474228 | REAC | Degradation of the extracellular matrix | 11056 | 748 | REAC:R-HSA-1474244 | REAC,REAC,REAC,REAC,REAC | COL6A1,TIMP2,COL6A2,COL5A1,FBN1 |
| query_1 | TRUE | 0.0316913 | 39 | 57 | 3 | 0.0526316 | 0.0769231 | REAC:R-HSA-140877 | REAC | Formation of Fibrin Clot (Clotting Cascade) | 11056 | 986 | REAC:R-HSA-109582 | REAC,REAC,REAC | F3,F2R,PROCR |
| query_1 | TRUE | 0.0322168 | 42 | 57 | 3 | 0.0526316 | 0.0714286 | REAC:R-HSA-419037 | REAC | NCAM1 interactions | 11056 | 1577 | REAC:R-HSA-375165 | REAC,REAC,REAC | COL6A1,COL6A2,COL5A1 |
| query_1 | TRUE | 0.0322168 | 516 | 57 | 9 | 0.1578947 | 0.0174419 | REAC:R-HSA-9006934 | REAC | Signaling by Receptor Tyrosine Kinases | 11056 | 2399 | REAC:R-HSA-162582 | REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC | CAV1,FOSL1,AXL,COL6A1,F3,TCIRG1,COL6A2,COL5A1,CYBA |
| query_1 | TRUE | 0.0338224 | 44 | 57 | 3 | 0.0526316 | 0.0681818 | REAC:R-HSA-8948216 | REAC | Collagen chain trimerization | 11056 | 401 | REAC:R-HSA-1650814 | REAC,REAC,REAC | COL6A1,COL6A2,COL5A1 |
| query_1 | TRUE | 0.0449470 | 675 | 57 | 10 | 0.1754386 | 0.0148148 | REAC:R-HSA-109582 | REAC | Hemostasis | 11056 | 1168 | REAC:0000000 | REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC,REAC | CAV1,SERPINE1,F3,PDE9A,ITGA5,EHD2,F2R,PROCR,L1CAM,SERPINA1 |
YOUR TURN: DISPLAY THE TOP 20 PATHWAYS FOR THE WIKIPATHWAY DATABASE (“WP”) and the down genes
enrichment_results <- gprofiler_results_DOWNgenes$result
##Get the name of the databases
unique(enrichment_results$source)
##Filter the dataframe to retrieve only the Reactome results
enrichment_results_WP <- enrichment_results %>%
XXXX(source == "XXXX")
##Now display the top 20 pathways
kable(head(enrichment_results_WP, n=20), caption = "Enrichment Result") %>%
kable_styling(bootstrap_options = c("striped", "condensed"), font_size = 10)Step 5 - Filter results by geneset size
Restrict the results to just the ones that have at maximum gene-set size of 1000 and a minimum gene-set size of 10 and with an adjusted pvalue of 0.05.-
- maximum gene-set size of 1000: will remove very generic pathway terms that are less informative like “regulation of gene expression”
- minimum gene-set size of 10: will remove small pathways where the pvalue is significant athough the overlap size is very small
- p_thres: aim is to retrieve only pathways that are significant under this adjusted pvalue threshold of 0.01 (1% false positive)
enrichment_results <- gprofiler_results_UPgenes$result
min_gs_size = 10
max_gs_size = 1000
p_thres = 0.01
# filter by params defined above
# by default we have set the max and min gs size to 250 and 3, respectively.
##the pvalue was stored as a character vector, we have to transform back into a numerical vector to be able to use it to filter values
enrichment_results$p_value <- as.numeric(enrichment_results$p_value)
enrichment_results_mxgssize_1000_min_10_adjp_0_01 <-
subset(enrichment_results,
term_size >= min_gs_size &
term_size <= max_gs_size &
p_value <= p_thres )
#
myrows = nrow(enrichment_results_mxgssize_1000_min_10_adjp_0_01)
print(glue("It results in {myrows} selected pathways with maximum gene-set size of {max_gs_size} and minimum gene-set of {min_gs_size} with an adjusted p-value of {p_thres} . "))## It results in 15 selected pathways with maximum gene-set size of 1000 and minimum gene-set of 10 with an adjusted p-value of 0.01 .YOUR TURN: using the enrichment results for the down genes, use gene-set size filters and FDR filters and display the results
In a second step, try different thresholds of maximum gene-set size and adjusted pvalue threshold: 500 and 250 for * max_gs_size and 0.05 or 0.001 for pvalue-
- Option 1
- max_gs_size = 500
- p_thres = 0.001
or
- Option 2
- max_gs_size = 250
- p_thres = 0.05
How many pathways do you get for each option (nrow())??
enrichment_results <- gprofiler_results_DOWNgenes$result
min_gs_size = XXXX
max_gs_size = XXXX
p_thres = XXXX
# filter by params defined above
# by default we have set the max and min gs size to 250 and 3, respectively.
##the pvalue was stored as a character vector, we have to transform back into a numerical vector to be able to use it to filter values
enrichment_results$p_value <- as.numeric(enrichment_results$p_value)
enrichment_results <- XXXXX
#
myrows = nrow(enrichment_results)
print(glue("It results in {myrows} selected pathways with maximum gene-set size of {max_gs_size} and minimum gene-set of {min_gs_size} with an adjusted p-value of {p_thres} . "))Step 6 - Run g:profiler with your own genesets (example using BaderLab genesets)
With regards to pathway sets there are two options when using g:Profiler -
- Use the genesets that are supplied by g:Profiler as we just did.
- Upload your own genesets.
The most common reasons for supplying your own genesets is the ability to use up to date annotations or in-house annotations that might not be available in the public sphere yet. It can be used to test the enrichment in a particular pathway of interest that might not be in the general pathway database but that you obtained for example by extracting the data from a published paper. You need to format this pathway in a .gmt file.
In this example, we will upload the Baderlab gene set file that contains multiple database sources. We filtered the gene-set prior to uploading the gene-set file to g:Profiler. The Baderlab gene set file combines different sources of databases and you can find more information at: https://baderlab.org/GeneSets
’‘’Why do we filter the geneset before running gprofiler’’’
The g:Profiler interface only allows for filtering genesets by size only after the analysis is complete. After the analysis is complete means the filtering is happening after Multiple hypothesis testing. Filtering prior to the analysis will generate more robust results because we exclude the uninformative large genesets prior to testing changing the sets that multiple hypothesis filtering will get rid of.
We have 1 example of gmt files with different filtering thresholds : * genesets less of equal than 500 genes
The format of the GMT file is described https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats#GMT:Gene_Matrix_Transposed_file_format.28.2A.gmt.29 and consists of rows with the following
- Name
- Description
- tab delimited list of genes a part of this geneset
Write out the gmt file with genenames
In order to use your own genesets with g:Profiler you need to upload the the file to g:Profiler server first. The function will return an ID that you need to specify in the organism parameter of the g:Profiler gost function call. It is done by using the upload_GMT_file() function which is included in the gprofiler2 package. Note: it took about 2 minutes to run this chunck of code.
dest_gmt_file_500 =file.path(data_dir, "data/Human_GOBP_AllPathways_filtered_500orLessV2.gmt")
custom_gmt_max500 <- upload_GMT_file(
gmtfile=dest_gmt_file_500)## Your custom annotations ID is gp__1EtK_ZBaC_KKk.
## You can use this ID as an 'organism' name in all the related enrichment tests against this custom source.## Just use: gost(my_genes, organism = 'gp__1EtK_ZBaC_KKk')UPgenes <- res$gene[res$logFC>0 & res$FDR <= 0.001]
DOWNgenes <- res$gene[res$logFC<0 & res$FDR <= 0.001]
gprofiler_results_custom_max500_UPgenes <- gost(query = UPgenes ,
significant=TRUE,
ordered_query = FALSE,
exclude_iea=TRUE,
correction_method = "fdr",
organism = custom_gmt_max500,
evcodes = TRUE
)## Detected custom GMT source requestenrichment_results_UPgenes <- gprofiler_results_custom_max500_UPgenes$result
enrichment_results_UPgenes$parents <- sapply(
enrichment_results_UPgenes$parents,
function(x) paste(x, collapse = ", ")
)
##Now display the top 20 pathways
kable(head(enrichment_results_UPgenes, n=20), caption = "Enrichment Result") %>%
kable_styling(bootstrap_options = c("striped", "condensed"), font_size = 10)| query | significant | p_value | term_size | query_size | intersection_size | precision | recall | term_id | source | term_name | effective_domain_size | source_order | parents | evidence_codes | intersection |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| query_1 | TRUE | 0.0000025 | 308 | 81 | 13 | 0.1604938 | 0.0422078 | EXTRACELLULAR MATRIX ORGANIZATION%REACTOME%R-HSA-1474244.5 | Human_GOBP_AllPathways_filtered_500orLessV2 | Extracellular matrix organization | 17710 | 3114 | X,X,X,X,X,X,X,X,X,X,X,X,X | COL6A1,SERPINE1,SNTB1,P3H2,TIMP2,P3H3,COL6A2,COL5A1,ITGA5,FBLN1,FBN1,LOX,BGN | |
| query_1 | TRUE | 0.0005125 | 66 | 81 | 6 | 0.0740741 | 0.0909091 | BETA1 INTEGRIN CELL SURFACE INTERACTIONS%PATHWAY INTERACTION DATABASE NCI-NATURE CURATED DATA%BETA1 INTEGRIN CELL SURFACE INTERACTIONS | Human_GOBP_AllPathways_filtered_500orLessV2 | Beta1 integrin cell surface interactions | 17710 | 1083 | X,X,X,X,X,X | COL6A1,COL6A2,TGFBI,COL5A1,ITGA5,FBN1 | |
| query_1 | TRUE | 0.0020101 | 89 | 81 | 6 | 0.0740741 | 0.0674157 | COLLAGEN FORMATION%REACTOME DATABASE ID RELEASE 97%1474290 | Human_GOBP_AllPathways_filtered_500orLessV2 | Collagen formation | 17710 | 2582 | X,X,X,X,X,X | COL6A1,P3H2,P3H3,COL6A2,COL5A1,LOX | |
| query_1 | TRUE | 0.0043121 | 62 | 81 | 5 | 0.0617284 | 0.0806452 | REGULATION OF COAGULATION%GOBP%GO:0050818 | Human_GOBP_AllPathways_filtered_500orLessV2 | regulation of coagulation | 17710 | 14756 | X,X,X,X,X | CAV1,SERPINE1,F3,F2R,PROCR | |
| query_1 | TRUE | 0.0050573 | 67 | 81 | 5 | 0.0617284 | 0.0746269 | COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES%REACTOME DATABASE ID RELEASE 97%1650814 | Human_GOBP_AllPathways_filtered_500orLessV2 | Collagen biosynthesis and modifying enzymes | 17710 | 3202 | X,X,X,X,X | COL6A1,P3H2,P3H3,COL6A2,COL5A1 | |
| query_1 | TRUE | 0.0067961 | 422 | 81 | 10 | 0.1234568 | 0.0236967 | TUBE DEVELOPMENT%GOBP%GO:0035295 | Human_GOBP_AllPathways_filtered_500orLessV2 | tube development | 17710 | 9957 | X,X,X,X,X,X,X,X,X,X | CAV1,SERPINE1,TGFBI,ITGA5,CCBE1,LOX,WNK4,ERRFI1,RAMP1,CLMP | |
| query_1 | TRUE | 0.0086573 | 40 | 81 | 4 | 0.0493827 | 0.1000000 | ELASTIC FIBRE FORMATION%REACTOME DATABASE ID RELEASE 97%1566948 | Human_GOBP_AllPathways_filtered_500orLessV2 | Elastic fibre formation | 17710 | 4843 | X,X,X,X | ITGA5,FBLN1,FBN1,LOX | |
| query_1 | TRUE | 0.0104444 | 283 | 81 | 8 | 0.0987654 | 0.0282686 | VASCULATURE DEVELOPMENT%GOBP%GO:0001944 | Human_GOBP_AllPathways_filtered_500orLessV2 | vasculature development | 17710 | 8945 | X,X,X,X,X,X,X,X | CAV1,SERPINE1,TGFBI,ITGA5,CCBE1,LOX,ERRFI1,RAMP1 | |
| query_1 | TRUE | 0.0137184 | 57 | 81 | 4 | 0.0493827 | 0.0701754 | COAGULATION PATHWAY%REACTOME DATABASE ID RELEASE 97%9769740 | Human_GOBP_AllPathways_filtered_500orLessV2 | Coagulation pathway | 17710 | 2450 | X,X,X,X | SERPINE1,F3,F2R,PROCR | |
| query_1 | TRUE | 0.0137184 | 98 | 81 | 5 | 0.0617284 | 0.0510204 | POSITIVE REGULATION OF MONOATOMIC ION TRANSMEMBRANE TRANSPORT%GOBP%GO:0034767 | Human_GOBP_AllPathways_filtered_500orLessV2 | positive regulation of monoatomic ion transmembrane transport | 17710 | 7284 | X,X,X,X,X | CAV1,F2R,FXYD5,TESC,FHL1 | |
| query_1 | TRUE | 0.0137184 | 50 | 81 | 4 | 0.0493827 | 0.0800000 | REGULATION OF CLOTTING CASCADE%REACTOME DATABASE ID RELEASE 97%9769739 | Human_GOBP_AllPathways_filtered_500orLessV2 | Regulation of clotting cascade | 17710 | 4167 | X,X,X,X | SERPINE1,F3,F2R,PROCR | |
| query_1 | TRUE | 0.0137184 | 57 | 81 | 4 | 0.0493827 | 0.0701754 | REGULATION OF BLOOD COAGULATION%GOBP%GO:0030193 | Human_GOBP_AllPathways_filtered_500orLessV2 | regulation of blood coagulation | 17710 | 17655 | X,X,X,X | CAV1,SERPINE1,F3,F2R | |
| query_1 | TRUE | 0.0137184 | 21 | 81 | 3 | 0.0370370 | 0.1428571 | POSITIVE REGULATION OF HEMOSTASIS%GOBP%GO:1900048 | Human_GOBP_AllPathways_filtered_500orLessV2 | positive regulation of hemostasis | 17710 | 15624 | X,X,X | SERPINE1,F3,F2R | |
| query_1 | TRUE | 0.0137184 | 240 | 81 | 7 | 0.0864198 | 0.0291667 | BLOOD VESSEL MORPHOGENESIS%GOBP%GO:0048514 | Human_GOBP_AllPathways_filtered_500orLessV2 | blood vessel morphogenesis | 17710 | 12918 | X,X,X,X,X,X,X | CAV1,SERPINE1,TGFBI,ITGA5,CCBE1,LOX,RAMP1 | |
| query_1 | TRUE | 0.0137184 | 94 | 81 | 5 | 0.0617284 | 0.0531915 | POSITIVE REGULATION OF CATION TRANSMEMBRANE TRANSPORT%GOBP%GO:1904064 | Human_GOBP_AllPathways_filtered_500orLessV2 | positive regulation of cation transmembrane transport | 17710 | 11849 | X,X,X,X,X | CAV1,F2R,FXYD5,TESC,FHL1 | |
| query_1 | TRUE | 0.0137184 | 56 | 81 | 4 | 0.0493827 | 0.0714286 | ASSEMBLY OF COLLAGEN FIBRILS AND OTHER MULTIMERIC STRUCTURES%REACTOME DATABASE ID RELEASE 97%2022090 | Human_GOBP_AllPathways_filtered_500orLessV2 | Assembly of collagen fibrils and other multimeric structures | 17710 | 4570 | X,X,X,X | COL6A1,COL6A2,COL5A1,LOX | |
| query_1 | TRUE | 0.0137184 | 21 | 81 | 3 | 0.0370370 | 0.1428571 | POSITIVE REGULATION OF BLOOD COAGULATION%GOBP%GO:0030194 | Human_GOBP_AllPathways_filtered_500orLessV2 | positive regulation of blood coagulation | 17710 | 18367 | X,X,X | SERPINE1,F3,F2R | |
| query_1 | TRUE | 0.0137184 | 52 | 81 | 4 | 0.0493827 | 0.0769231 | COLLAGEN METABOLIC PROCESS%GOBP%GO:0032963 | Human_GOBP_AllPathways_filtered_500orLessV2 | collagen metabolic process | 17710 | 16593 | X,X,X,X | P3H2,P3H3,COL5A1,MRC2 | |
| query_1 | TRUE | 0.0139108 | 58 | 81 | 4 | 0.0493827 | 0.0689655 | REGULATION OF HEMOSTASIS%GOBP%GO:1900046 | Human_GOBP_AllPathways_filtered_500orLessV2 | regulation of hemostasis | 17710 | 7751 | X,X,X,X | CAV1,SERPINE1,F3,F2R | |
| query_1 | TRUE | 0.0141423 | 23 | 81 | 3 | 0.0370370 | 0.1304348 | POSITIVE REGULATION OF COAGULATION%GOBP%GO:0050820 | Human_GOBP_AllPathways_filtered_500orLessV2 | positive regulation of coagulation | 17710 | 12228 | X,X,X | SERPINE1,F3,F2R |
YOUR TURN: perform the same task with the down genes
gprofiler_results_custom_max500_DNgenes <- gost(
query = XXXX ,
significant=XXXX,
ordered_query = XXXX,
exclude_iea=XXXX,
correction_method = "XXX",
organism = XXXX,
evcodes = XXXX
)
enrichment_results_DNgenes <- gprofiler_results_custom_max500_DNgenes$result
enrichment_results_DNgenes$parents <- sapply(
enrichment_results_DNgenes$parents,
function(x) paste(x, collapse = ", ")
)
kable(head(enrichment_results_DNgenes, n=20), caption = "Enrichment Result") %>%
kable_styling(bootstrap_options = c("striped", "condensed"), font_size = 10)
# write.csv( enrichment_results_DNgenes, "results/enrichment_results_gprofiler_DNgenes.csv")Exercise 2 - visualization
We will use clusterProfiler style visualization. ### Step 1 - Bar plot
##Calculate a score (small very significant adjusted p value will turn into a large score)
enrichment_results_UPgenes <- gprofiler_results_custom_max500_UPgenes$result
enrichment_results_UPgenes$score = -log10(as.numeric(enrichment_results_UPgenes$p_value))
enrichment <- enrichment_results_UPgenes %>%
as_tibble() %>%
arrange(desc(score)) %>%
slice_head(n = 20)
p = ggplot(enrichment, aes(reorder(term_name, score), score)) +
geom_col(aes(fill = score), width=0.8) +
scale_fill_gradient(low = "#C1F6F8", high = "#00BFC4")+
coord_flip() +
labs(x="", y="Normalized Enrichment Score",title="Pathway enrichment analysis ") +
theme_minimal()
p
Step 2 - Dot plot with segment
enrichmentTop20 <- enrichment %>%
arrange(desc(score)) %>%
slice_head(n = 20)
p = ggplot(enrichmentTop20, aes(
y = fct_reorder(term_name, score),
x = score,
color = score,
size = score
)) +
geom_point() +
scale_color_gradientn(
colours = c("#56A443", "#8AA443", "#f7ca64"),
trans = "log10",
guide = guide_colorbar(reverse = FALSE, order = 1)
) +
scale_size_continuous(range = c(4, 8)) +
theme_bw(base_size = 12) +
xlab("score -log10(adj pvalue)") +
ylab(NULL) +
ggtitle("Pathway Enrichment Analysis") +
theme(
axis.text.y = element_text(size = 10)
)
p
Step 4 - Hierarchical Cluster
The problem with the earlier visualization, which represents the top 20 pathways, is that it does not account for the overlap between gene sets. As you can observe, many pathways have similar names because they share a substantial number of genes. To provide a more accurate and efficient representation of the results, and to obtain a quicker and clearer overview, it would be more appropriate to aggregate and summarize pathways with substantial gene overlap, thereby reducing redundancy.
Our first approach is to use hierarchical clustering based on a distance matrix that measures the similarity between pathways according to the number of genes they have in common, allowing us to visualize their relationships.
enrichment_results_UPgenes <- gprofiler_results_custom_max500_UPgenes$result
enrichment_results_UPgenes$score = -log10(as.numeric(enrichment_results_UPgenes$p_value))
# Top 20 enriched terms
enrichmentTop20 <- enrichment_results_UPgenes %>%
arrange(desc(score)) %>%
slice_head(n = 20)Calculating the distance matrix using the jaccard coefficient
- Many pathways share genes, so similar biological processes can appear under different pathway names. To group similar pathways, we look at how many genes they have in common.
- The Jaccard coefficient measures the shared genes relative to all unique genes in the two pathways, from 0 (no overlap) to 1 (identical).
- The overlap coefficient measures the shared genes relative to the smaller pathway and shows when one pathway is mostly contained within another, from 0 (no overlap) to 1 (smaller pathway contained into larger one).
Figure 3: Jaccard and overlap coefficients for two gene sets.
We convert these scores into distances (1 − score) and use them to group similar pathways into broader biological themes.
##Illustration with Jaccard Coefficient
# Convert genes to lists
# (the intersection column, returned by g:Profiler with evcodes = TRUE, contains
# the genes of our list found in each pathway, separated by a comma)
gene_lists <- strsplit(enrichment_results_UPgenes$intersection, ",")
# Calculate Jaccard similarity matrix
n <- length(gene_lists)
jaccard_mat <- matrix(0, n, n)
rownames(jaccard_mat) <- enrichment_results_UPgenes$term_name
colnames(jaccard_mat) <- enrichment_results_UPgenes$term_name
for (i in seq_len(n)) {
for (j in seq_len(n)) {
intersection <- length(intersect(gene_lists[[i]], gene_lists[[j]]))
union <- length(union(gene_lists[[i]], gene_lists[[j]]))
jaccard_mat[i, j] <- intersection / union
}
}
labs <- enrichment_results_UPgenes$term_name
labs <- ifelse(nchar(labs) > 30, paste0(substr(labs, 1, 27), "..."), labs)
rownames(jaccard_mat) <- colnames(jaccard_mat) <- make.unique(labs)
keep <- sample(seq_len(nrow(jaccard_mat)), min(20, nrow(jaccard_mat)))
jaccard_sub <- jaccard_mat[keep, keep]
# Heatmap
p <- pheatmap(
jaccard_sub,
color = colorRampPalette(c("white", "#C6DBEF", "#6BAED6", "#2171B5", "#08306B"))(100),
breaks = seq(0, 1, length.out = 101),
display_numbers = F,
#number_format = "%.2f",
#number_color = "black",
#fontsize_number = 8,
clustering_distance_rows = as.dist(1 - jaccard_sub),
clustering_distance_cols = as.dist(1 - jaccard_sub),
border_color = "grey90",
fontsize_row = 7,
fontsize_col = 7,
angle_col = 45,
main = "Jaccard coefficient (20 random gene sets)"
)
##Illustration with Overlap Coefficient
# Calculate Overlap coefficient matrix (reuses gene_lists and n from above)
overlap_mat <- matrix(0, n, n)
for (i in seq_len(n)) {
for (j in seq_len(n)) {
intersection <- length(intersect(gene_lists[[i]], gene_lists[[j]]))
smaller <- min(length(gene_lists[[i]]), length(gene_lists[[j]]))
overlap_mat[i, j] <- intersection / smaller
}
}
rownames(overlap_mat) <- colnames(overlap_mat) <- make.unique(labs)
# Use the same 20 gene sets as the Jaccard heatmap so the two can be compared
overlap_sub <- overlap_mat[keep, keep]
# Heatmap
p <- pheatmap(
overlap_sub,
color = colorRampPalette(c("white", "#C6DBEF", "#6BAED6", "#2171B5", "#08306B"))(100),
breaks = seq(0, 1, length.out = 101),
display_numbers = F,
clustering_distance_rows = as.dist(1 - overlap_sub),
clustering_distance_cols = as.dist(1 - overlap_sub),
border_color = "grey90",
fontsize_row = 7,
fontsize_col = 7,
angle_col = 45,
main = "Overlap coefficient (20 random gene sets)"
)
##Hierarchical Clustering
plot_hclust_enrichment: helper function to cluster the top enriched pathways by gene overlap (Jaccard distance) and plot them next to the dendrogram
plot_hclust_enrichment = function(enrichment_results, n_top = 20, n_clusters = 8) {
enrichment_results$score <- -log10(as.numeric(enrichment_results$p_value))
# Top n_top enriched terms
enrichmentTop <- enrichment_results %>%
arrange(desc(score)) %>%
slice_head(n = n_top)
top_terms <- enrichmentTop %>%
arrange(as.numeric(p_value))
# Jaccard similarity matrix
# (the intersection column contains the genes of our list found in each pathway)
gene_lists <- strsplit(top_terms$intersection, ",")
n <- length(gene_lists)
jaccard_mat <- matrix(0, n, n,
dimnames = list(top_terms$term_name, top_terms$term_name))
for (i in seq_len(n)) {
for (j in seq_len(n)) {
inter <- length(intersect(gene_lists[[i]], gene_lists[[j]]))
uni <- length(union(gene_lists[[i]], gene_lists[[j]]))
jaccard_mat[i, j] <- inter / uni
}
}
# Hierarchical clustering
dist_mat <- as.dist(1 - jaccard_mat)
hc <- hclust(dist_mat, method = "average")
term_order <- hc$labels[hc$order]
# The first level is drawn at the bottom of the bar plot, matching the
# first dendrogram leaf (y = 1)
top_terms$term_name <- factor(top_terms$term_name, levels = term_order)
top_terms$neglog10p <- -log10(as.numeric(top_terms$p_value))
# Cut tree into k branches (cannot have more branches than terms)
k <- min(n_clusters, n)
clusters <- cutree(hc, k = k)
top_terms$cluster <- factor(clusters[as.character(top_terms$term_name)])
# Bar plot (full, coloured by cluster)
base_colours <- c("#1D9E75", "#3A7DC9", "#E07B39", "#9B59B6",
"#E74C3C", "#F1C40F", "#1ABC9C", "#E91E63")
# Interpolate extra colours when more than 8 clusters are requested
if (k <= length(base_colours)) {
cluster_colours <- base_colours
} else {
cluster_colours <- colorRampPalette(base_colours)(k)
}
p_bar <- ggplot(top_terms, aes(x = neglog10p, y = term_name, fill = cluster)) +
geom_col(alpha = 0.85, width = 0.7) +
scale_y_discrete(expand = c(0, 0.5)) +
scale_fill_manual(values = cluster_colours[seq_len(k)],
name = "Branch") +
labs(x = expression(-log[10](p-value)), y = NULL) +
theme_bw(base_size = 10) +
theme(
panel.grid.major.y = element_blank(),
axis.text.y = element_text(size = 8),
legend.position = "right"
)
# Dendrogram
dend_data <- dendro_data(hc, type = "rectangle")
seg <- dend_data$segments
p_dend <- ggplot(seg) +
geom_segment(aes(x = y, y = x, xend = yend, yend = xend),
colour = "grey40", linewidth = 0.45) +
scale_x_reverse(expand = expansion(mult = c(0.02, 0))) +
scale_y_continuous(limits = c(0.5, n + 0.5), expand = c(0, 0)) +
labs(x = "Jaccard distance") +
theme_void() +
theme(
axis.title.x = element_text(size = 8, colour = "grey50", margin = margin(t = 4)),
axis.text.x = element_text(size = 7, colour = "grey60"),
plot.margin = margin(0, 4, 0, 0)
)
# Combined plot
p_combined <- p_dend + p_bar +
plot_layout(widths = c(3, 1))
# Keep only the top hit per branch
top_per_cluster <- top_terms %>%
group_by(cluster) %>%
slice_min(as.numeric(p_value), n = 1, with_ties = FALSE) %>%
ungroup() %>%
arrange(as.numeric(p_value)) %>%
mutate(term_name = factor(as.character(term_name),
levels = rev(as.character(term_name))))
# Simple barplot — one representative term per branch
p_simple <- ggplot(top_per_cluster,
aes(x = neglog10p, y = term_name, fill = cluster)) +
geom_col(alpha = 0.88, width = 0.6) +
geom_text(aes(label = sprintf("p = %.2e", as.numeric(p_value))),
hjust = -0.1, size = 3, colour = "grey30") +
scale_fill_manual(values = cluster_colours[seq_len(k)],
name = "Branch") +
scale_x_continuous(expand = expansion(mult = c(0, 0.25))) +
labs(
x = expression(-log[10](p-value)),
y = NULL,
title = "Top term per branch"
) +
theme_bw(base_size = 11) +
theme(
panel.grid.major.y = element_blank(),
axis.text.y = element_text(size = 9),
legend.position = "bottom",
plot.title = element_text(size = 11, face = "bold")
)
# n_terms and n_clusters are returned so the figure height can be set
# from them in the chunk options
return(list(combined = p_combined, top_per_branch = p_simple,
n_terms = n, n_clusters = k))
}Running the function on the top 20 pathways
enrichment_results_UPgenes <- gprofiler_results_custom_max500_UPgenes$result
hc_plots <- plot_hclust_enrichment(enrichment_results = enrichment_results_UPgenes, n_top = 20, n_clusters = 10)The figure height grows with the number of pathways (combined plot) and the number of clusters (top term per branch), so the bars do not overlap.
