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
)
## Loaded Tensorflow version 2.7.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.6765993
g2_term_3 -0.6015774
g2_term_2 -0.5822843
g2_term_4 -0.4234689
g2_term_8 -0.4132009
g1_term_2 0.3973174

…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")
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

treeplot(obj) +
  scale_colour_continuous(name = "Enrichment Score")
## 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.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server x64 (build 20348)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] pareg_1.2.0           tfprobability_0.15.1  tensorflow_2.9.0     
##  [4] enrichplot_1.18.0     ComplexHeatmap_2.14.0 forcats_0.5.2        
##  [7] stringr_1.4.1         dplyr_1.0.10          purrr_0.3.5          
## [10] readr_2.1.3           tidyr_1.2.1           tibble_3.1.8         
## [13] tidyverse_1.3.2       ggraph_2.1.0          ggplot2_3.4.0        
## [16] BiocStyle_2.26.0     
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2             reticulate_1.26        tidyselect_1.2.0      
##   [4] RSQLite_2.2.18         AnnotationDbi_1.60.0   BiocParallel_1.32.1   
##   [7] scatterpie_0.1.8       munsell_0.5.0          codetools_0.2-18      
##  [10] future_1.29.0          withr_2.5.0            keras_2.9.0           
##  [13] colorspace_2.0-3       GOSemSim_2.24.0        Biobase_2.58.0        
##  [16] highr_0.9              knitr_1.40             stats4_4.2.2          
##  [19] DOSE_3.24.1            listenv_0.8.0          labeling_0.4.2        
##  [22] GenomeInfoDbData_1.2.9 matrixLaplacian_1.0    polyclip_1.10-4       
##  [25] bit64_4.0.5            farver_2.1.1           parallelly_1.32.1     
##  [28] vctrs_0.5.0            treeio_1.22.0          generics_0.1.3        
##  [31] xfun_0.34              timechange_0.1.1       R6_2.5.1              
##  [34] doParallel_1.0.17      GenomeInfoDb_1.34.3    clue_0.3-62           
##  [37] graphlayouts_0.8.3     bitops_1.0-7           cachem_1.0.6          
##  [40] fgsea_1.24.0           gridGraphics_0.5-1     assertthat_0.2.1      
##  [43] scales_1.2.1           googlesheets4_1.0.1    gtable_0.3.1          
##  [46] Cairo_1.6-0            globals_0.16.1         tidygraph_1.2.2       
##  [49] rlang_1.0.6            zeallot_0.1.0          scatterplot3d_0.3-42  
##  [52] GlobalOptions_0.1.2    splines_4.2.2          lazyeval_0.2.2        
##  [55] gargle_1.2.1           broom_1.0.1            BiocManager_1.30.19   
##  [58] yaml_2.3.6             reshape2_1.4.4         modelr_0.1.10         
##  [61] backports_1.4.1        qvalue_2.30.0          tools_4.2.2           
##  [64] bookdown_0.30          ggplotify_0.1.0        ellipsis_0.3.2        
##  [67] jquerylib_0.1.4        RColorBrewer_1.1-3     proxy_0.4-27          
##  [70] BiocGenerics_0.44.0    Rcpp_1.0.9             plyr_1.8.8            
##  [73] progress_1.2.2         base64enc_0.1-3        zlibbioc_1.44.0       
##  [76] RCurl_1.98-1.9         prettyunits_1.1.1      GetoptLong_1.0.5      
##  [79] viridis_0.6.2          cowplot_1.1.1          S4Vectors_0.36.0      
##  [82] haven_2.5.1            ggrepel_0.9.2          cluster_2.1.4         
##  [85] fs_1.5.2               furrr_0.3.1            magrittr_2.0.3        
##  [88] magick_2.7.3           data.table_1.14.4      circlize_0.4.15       
##  [91] reprex_2.0.2           googledrive_2.0.0      whisker_0.4           
##  [94] ggnewscale_0.4.8       matrixStats_0.62.0     hms_1.1.2             
##  [97] patchwork_1.1.2        evaluate_0.18          HDO.db_0.99.1         
## [100] readxl_1.4.1           IRanges_2.32.0         gridExtra_2.3         
## [103] shape_1.4.6            tfruns_1.5.1           compiler_4.2.2        
## [106] crayon_1.5.2           shadowtext_0.1.2       htmltools_0.5.3       
## [109] ggfun_0.0.8            tzdb_0.3.0             aplot_0.1.8           
## [112] lubridate_1.9.0        DBI_1.1.3              tweenr_2.0.2          
## [115] dbplyr_2.2.1           MASS_7.3-58.1          Matrix_1.5-3          
## [118] cli_3.4.1              parallel_4.2.2         igraph_1.3.5          
## [121] pkgconfig_2.0.3        xml2_1.3.3             foreach_1.5.2         
## [124] ggtree_3.6.2           bslib_0.4.1            XVector_0.38.0        
## [127] rvest_1.0.3            yulab.utils_0.0.5      digest_0.6.30         
## [130] Biostrings_2.66.0      rmarkdown_2.18         cellranger_1.1.0      
## [133] fastmatch_1.1-3        tidytree_0.4.1         nloptr_2.0.3          
## [136] rjson_0.2.21           lifecycle_1.0.3        nlme_3.1-160          
## [139] jsonlite_1.8.3         viridisLite_0.4.1      fansi_1.0.3           
## [142] pillar_1.8.1           lattice_0.20-45        KEGGREST_1.38.0       
## [145] fastmap_1.1.0          httr_1.4.4             GO.db_3.16.0          
## [148] glue_1.6.2             png_0.1-7              iterators_1.0.14      
## [151] bit_4.0.4              ggforce_0.4.1          stringi_1.7.8         
## [154] sass_0.4.2             blob_1.2.3             memoise_2.0.1         
## [157] ape_5.6-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.