Contents

1 Overview

This vignette is an introduction to the usage of pareg. It estimates pathway enrichment scores by regressing differential expression p-values of all genes considered in an experiment on their membership to a set of biological pathways. These scores are computed using a regularized generalized linear model with LASSO and network regularization terms. The network regularization term is based on a pathway similarity matrix (e.g., defined by Jaccard similarity) and thus classifies this method as a modular enrichment analysis tool (Huang, Sherman, and Lempicki 2009).

2 Installation

if (!require("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("pareg")

3 Load required packages

We start our analysis by loading the pareg package and other required libraries.

library(ggraph)
library(tidyverse)
library(ComplexHeatmap)
library(enrichplot)

library(pareg)

set.seed(42)

4 Introductory example

4.1 Generate pathway database

For the sake of this introductory example, we generate a synthetic pathway database with a pronounced clustering of pathways.

group_num <- 2
pathways_from_group <- 10

gene_groups <- purrr::map(seq(1, group_num), function(group_idx) {
  glue::glue("g{group_idx}_gene_{seq_len(15)}")
})
genes_bg <- paste0("bg_gene_", seq(1, 50))

df_terms <- purrr::imap_dfr(
  gene_groups,
  function(current_gene_list, gene_list_idx) {
    purrr::map_dfr(seq_len(pathways_from_group), function(pathway_idx) {
      data.frame(
        term = paste0("g", gene_list_idx, "_term_", pathway_idx),
        gene = c(
          sample(current_gene_list, 10, replace = FALSE),
          sample(genes_bg, 10, replace = FALSE)
        )
      )
    })
  }
)

df_terms %>%
  sample_n(5)
##        term       gene
## 1 g1_term_9 g1_gene_12
## 2 g1_term_5  g1_gene_7
## 3 g2_term_2  g2_gene_2
## 4 g1_term_3 bg_gene_47
## 5 g1_term_8  g1_gene_1

4.2 Term similarities

Before starting the actual enrichment estimation, we compute pairwise pathway similarities with pareg’s helper function.

mat_similarities <- compute_term_similarities(
  df_terms,
  similarity_function = jaccard
)

hist(mat_similarities, xlab = "Term similarity")

We can see a clear clustering of pathways.

Heatmap(
  mat_similarities,
  name = "Similarity",
  col = circlize::colorRamp2(c(0, 1), c("white", "black"))
)

4.3 Create synthetic study

We then select a subset of pathways to be activated. In a performance evaluation, these would be considered to be true positives.

active_terms <- similarity_sample(mat_similarities, 5)
active_terms
## [1] "g2_term_6" "g2_term_3" "g2_term_3" "g2_term_2" "g2_term_8"

The genes contained in the union of active pathways are considered to be differentially expressed.

de_genes <- df_terms %>%
  filter(term %in% active_terms) %>%
  distinct(gene) %>%
  pull(gene)

other_genes <- df_terms %>%
  distinct(gene) %>%
  pull(gene) %>%
  setdiff(de_genes)

The p-values of genes considered to be differentially expressed are sampled from a Beta distribution centered at \(0\). The p-values for all other genes are drawn from a Uniform distribution.

df_study <- data.frame(
  gene = c(de_genes, other_genes),
  pvalue = c(rbeta(length(de_genes), 0.1, 1), rbeta(length(other_genes), 1, 1)),
  in_study = c(
    rep(TRUE, length(de_genes)),
    rep(FALSE, length(other_genes))
  )
)

table(
  df_study$pvalue <= 0.05,
  df_study$in_study, dnn = c("sig. p-value", "in study")
)
##             in study
## sig. p-value FALSE TRUE
##        FALSE    34   17
##        TRUE      1   28

4.4 Enrichment analysis

Finally, we compute pathway enrichment scores.

fit <- pareg(
  df_study %>% select(gene, pvalue),
  df_terms,
  network_param = 1, term_network = mat_similarities
)
## + /home/biocbuild/.cache/R/basilisk/1.14.0/0/bin/conda 'create' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.14.0/pareg/1.6.0/pareg' 'python=3.9.12' '--quiet' '-c' 'anaconda'
## + /home/biocbuild/.cache/R/basilisk/1.14.0/0/bin/conda 'install' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.14.0/pareg/1.6.0/pareg' 'python=3.9.12'
## + /home/biocbuild/.cache/R/basilisk/1.14.0/0/bin/conda 'install' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.14.0/pareg/1.6.0/pareg' '-c' 'anaconda' 'python=3.9.12' 'tensorflow=2.10.0' 'tensorflow-probability=0.14.0'

The results can be exported to a dataframe for further processing…

fit %>%
  as.data.frame() %>%
  arrange(desc(abs(enrichment))) %>%
  head() %>%
  knitr::kable()
term enrichment
g2_term_6 -0.6759553
g2_term_3 -0.6003887
g2_term_2 -0.5817953
g2_term_4 -0.4232832
g2_term_8 -0.4122331
g1_term_2 0.3979995

…and also visualized in a pathway network view.

plot(fit, min_similarity = 0.1)

To provide a wider range of visualization options, the result can be transformed into an object which is understood by the functions of the enrichplot package.

obj <- as_enrichplot_object(fit)

dotplot(obj) +
  scale_colour_continuous(name = "Enrichment Score")

treeplot(obj) +
  scale_colour_continuous(name = "Enrichment Score")
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

4.5 Session information

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
## LAPACK: /home/biocbuild/.cache/R/basilisk/1.14.0/pareg/1.6.0/pareg/lib/libmkl_rt.so.1;  LAPACK version 3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] pareg_1.6.0           tfprobability_0.15.1  tensorflow_2.14.0    
##  [4] enrichplot_1.22.0     ComplexHeatmap_2.18.0 lubridate_1.9.3      
##  [7] forcats_1.0.0         stringr_1.5.0         dplyr_1.1.3          
## [10] purrr_1.0.2           readr_2.1.4           tidyr_1.3.0          
## [13] tibble_3.2.1          tidyverse_2.0.0       ggraph_2.1.0         
## [16] ggplot2_3.4.4         BiocStyle_2.30.0     
## 
## loaded via a namespace (and not attached):
##   [1] fs_1.6.3                      matrixStats_1.0.0            
##   [3] bitops_1.0-7                  devtools_2.4.5               
##   [5] HDO.db_0.99.1                 httr_1.4.7                   
##   [7] RColorBrewer_1.1-3            doParallel_1.0.17            
##   [9] profvis_0.3.8                 tools_4.3.1                  
##  [11] doRNG_1.8.6                   utf8_1.2.4                   
##  [13] R6_2.5.1                      lazyeval_0.2.2               
##  [15] GetoptLong_1.0.5              urlchecker_1.0.1             
##  [17] withr_2.5.1                   prettyunits_1.2.0            
##  [19] gridExtra_2.3                 cli_3.6.1                    
##  [21] Biobase_2.62.0                Cairo_1.6-1                  
##  [23] scatterpie_0.2.1              labeling_0.4.3               
##  [25] sass_0.4.7                    proxy_0.4-27                 
##  [27] yulab.utils_0.1.0             DOSE_3.28.0                  
##  [29] parallelly_1.36.0             sessioninfo_1.2.2            
##  [31] RSQLite_2.3.1                 generics_0.1.3               
##  [33] gridGraphics_0.5-1            shape_1.4.6                  
##  [35] GO.db_3.18.0                  Matrix_1.6-1.1               
##  [37] fansi_1.0.5                   S4Vectors_0.40.0             
##  [39] logger_0.2.2                  lifecycle_1.0.3              
##  [41] whisker_0.4.1                 yaml_2.3.7                   
##  [43] qvalue_2.34.0                 BiocFileCache_2.10.0         
##  [45] blob_1.2.4                    promises_1.2.1               
##  [47] crayon_1.5.2                  dir.expiry_1.10.0            
##  [49] miniUI_0.1.1.1                lattice_0.22-5               
##  [51] cowplot_1.1.1                 KEGGREST_1.42.0              
##  [53] magick_2.8.1                  zeallot_0.1.0                
##  [55] pillar_1.9.0                  knitr_1.44                   
##  [57] fgsea_1.28.0                  rjson_0.2.21                 
##  [59] future.apply_1.11.0           codetools_0.2-19             
##  [61] fastmatch_1.1-4               glue_1.6.2                   
##  [63] doFuture_1.0.0                ggfun_0.1.3                  
##  [65] data.table_1.14.8             remotes_2.4.2.1              
##  [67] vctrs_0.6.4                   png_0.1-8                    
##  [69] treeio_1.26.0                 gtable_0.3.4                 
##  [71] cachem_1.0.8                  xfun_0.40                    
##  [73] mime_0.12                     tidygraph_1.2.3              
##  [75] iterators_1.0.14              interactiveDisplayBase_1.40.0
##  [77] ellipsis_0.3.2                nlme_3.1-163                 
##  [79] ggtree_3.10.0                 usethis_2.2.2                
##  [81] bit64_4.0.5                   progress_1.2.2               
##  [83] filelock_1.0.2                GenomeInfoDb_1.38.0          
##  [85] bslib_0.5.1                   colorspace_2.1-0             
##  [87] BiocGenerics_0.48.0           DBI_1.1.3                    
##  [89] tidyselect_1.2.0              processx_3.8.2               
##  [91] bit_4.0.5                     compiler_4.3.1               
##  [93] curl_5.1.0                    basilisk.utils_1.14.0        
##  [95] bookdown_0.36                 shadowtext_0.1.2             
##  [97] scales_1.2.1                  keras_2.13.0                 
##  [99] callr_3.7.3                   tfruns_1.5.1                 
## [101] rappdirs_0.3.3                digest_0.6.33                
## [103] rmarkdown_2.25                basilisk_1.14.0              
## [105] XVector_0.42.0                htmltools_0.5.6.1            
## [107] pkgconfig_2.0.3               base64enc_0.1-3              
## [109] dbplyr_2.3.4                  fastmap_1.1.1                
## [111] rlang_1.1.1                   GlobalOptions_0.1.2          
## [113] htmlwidgets_1.6.2             shiny_1.7.5.1                
## [115] farver_2.1.1                  jquerylib_0.1.4              
## [117] jsonlite_1.8.7                BiocParallel_1.36.0          
## [119] GOSemSim_2.28.0               RCurl_1.98-1.12              
## [121] magrittr_2.0.3                GenomeInfoDbData_1.2.11      
## [123] ggplotify_0.1.2               patchwork_1.1.3              
## [125] munsell_0.5.0                 Rcpp_1.0.11                  
## [127] ggnewscale_0.4.9              ape_5.7-1                    
## [129] viridis_0.6.4                 reticulate_1.34.0            
## [131] stringi_1.7.12                zlibbioc_1.48.0              
## [133] MASS_7.3-60                   AnnotationHub_3.10.0         
## [135] plyr_1.8.9                    pkgbuild_1.4.2               
## [137] parallel_4.3.1                HPO.db_0.99.2                
## [139] listenv_0.9.0                 ggrepel_0.9.4                
## [141] Biostrings_2.70.0             graphlayouts_1.0.1           
## [143] splines_4.3.1                 hms_1.1.3                    
## [145] circlize_0.4.15               ps_1.7.5                     
## [147] igraph_1.5.1                  rngtools_1.5.2               
## [149] reshape2_1.4.4                stats4_4.3.1                 
## [151] pkgload_1.3.3                 BiocVersion_3.18.0           
## [153] evaluate_0.22                 BiocManager_1.30.22          
## [155] nloptr_2.0.3                  tzdb_0.4.0                   
## [157] foreach_1.5.2                 tweenr_2.0.2                 
## [159] httpuv_1.6.12                 polyclip_1.10-6              
## [161] future_1.33.0                 clue_0.3-65                  
## [163] ggforce_0.4.1                 xtable_1.8-4                 
## [165] tidytree_0.4.5                MPO.db_0.99.7                
## [167] later_1.3.1                   viridisLite_0.4.2            
## [169] aplot_0.2.2                   memoise_2.0.1                
## [171] AnnotationDbi_1.64.0          IRanges_2.36.0               
## [173] cluster_2.1.4                 timechange_0.2.0             
## [175] globals_0.16.2

References

Huang, Da Wei, Brad T Sherman, and Richard A Lempicki. 2009. “Bioinformatics Enrichment Tools: Paths Toward the Comprehensive Functional Analysis of Large Gene Lists.” Nucleic Acids Research 37 (1): 1–13.