# Get started

Kim Philipp Jablonski

#### 26 April 2022

#### Package

pareg 1.0.0

# 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.0 RC (2022-04-19 r82224 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.0.0           tfprobability_0.15.0  tensorflow_2.8.0
    ##  [4] enrichplot_1.16.0     ComplexHeatmap_2.12.0 forcats_0.5.1
    ##  [7] stringr_1.4.0         dplyr_1.0.8           purrr_0.3.4
    ## [10] readr_2.1.2           tidyr_1.2.0           tibble_3.1.6
    ## [13] tidyverse_1.3.1       ggraph_2.0.5          ggplot2_3.3.5
    ## [16] BiocStyle_2.24.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] utf8_1.2.2             reticulate_1.24        tidyselect_1.1.2
    ##   [4] RSQLite_2.2.12         AnnotationDbi_1.58.0   BiocParallel_1.30.0
    ##   [7] scatterpie_0.1.7       munsell_0.5.0          codetools_0.2-18
    ##  [10] future_1.25.0          withr_2.5.0            keras_2.8.0
    ##  [13] colorspace_2.0-3       GOSemSim_2.22.0        Biobase_2.56.0
    ##  [16] highr_0.9              knitr_1.38             rstudioapi_0.13
    ##  [19] stats4_4.2.0           DOSE_3.22.0            listenv_0.8.0
    ##  [22] labeling_0.4.2         GenomeInfoDbData_1.2.8 matrixLaplacian_1.0
    ##  [25] polyclip_1.10-0        bit64_4.0.5            farver_2.1.0
    ##  [28] rprojroot_2.0.3        parallelly_1.31.1      vctrs_0.4.1
    ##  [31] treeio_1.20.0          generics_0.1.2         xfun_0.30
    ##  [34] R6_2.5.1               doParallel_1.0.17      GenomeInfoDb_1.32.0
    ##  [37] clue_0.3-60            graphlayouts_0.8.0     bitops_1.0-7
    ##  [40] cachem_1.0.6           fgsea_1.22.0           gridGraphics_0.5-1
    ##  [43] assertthat_0.2.1       scales_1.2.0           gtable_0.3.0
    ##  [46] Cairo_1.5-15           globals_0.14.0         tidygraph_1.2.1
    ##  [49] rlang_1.0.2            zeallot_0.1.0          scatterplot3d_0.3-41
    ##  [52] GlobalOptions_0.1.2    splines_4.2.0          lazyeval_0.2.2
    ##  [55] broom_0.8.0            BiocManager_1.30.17    yaml_2.3.5
    ##  [58] reshape2_1.4.4         modelr_0.1.8           backports_1.4.1
    ##  [61] qvalue_2.28.0          tools_4.2.0            bookdown_0.26
    ##  [64] ggplotify_0.1.0        ellipsis_0.3.2         jquerylib_0.1.4
    ##  [67] RColorBrewer_1.1-3     proxy_0.4-26           BiocGenerics_0.42.0
    ##  [70] Rcpp_1.0.8.3           plyr_1.8.7             base64enc_0.1-3
    ##  [73] progress_1.2.2         zlibbioc_1.42.0        RCurl_1.98-1.6
    ##  [76] prettyunits_1.1.1      GetoptLong_1.0.5       viridis_0.6.2
    ##  [79] S4Vectors_0.34.0       haven_2.5.0            ggrepel_0.9.1
    ##  [82] cluster_2.1.3          here_1.0.1             fs_1.5.2
    ##  [85] furrr_0.2.3            magrittr_2.0.3         magick_2.7.3
    ##  [88] data.table_1.14.2      DO.db_2.9              circlize_0.4.14
    ##  [91] reprex_2.0.1           whisker_0.4            ggnewscale_0.4.7
    ##  [94] matrixStats_0.62.0     hms_1.1.1              patchwork_1.1.1
    ##  [97] evaluate_0.15          readxl_1.4.0           IRanges_2.30.0
    ## [100] gridExtra_2.3          shape_1.4.6            tfruns_1.5.0
    ## [103] compiler_4.2.0         crayon_1.5.1           shadowtext_0.1.2
    ## [106] htmltools_0.5.2        ggfun_0.0.6            tzdb_0.3.0
    ## [109] aplot_0.1.3            lubridate_1.8.0        DBI_1.1.2
    ## [112] tweenr_1.0.2           dbplyr_2.1.1           rappdirs_0.3.3
    ## [115] MASS_7.3-57            Matrix_1.4-1           cli_3.3.0
    ## [118] parallel_4.2.0         igraph_1.3.1           pkgconfig_2.0.3
    ## [121] xml2_1.3.3             foreach_1.5.2          ggtree_3.4.0
    ## [124] bslib_0.3.1            XVector_0.36.0         rvest_1.0.2
    ## [127] yulab.utils_0.0.4      digest_0.6.29          Biostrings_2.64.0
    ## [130] rmarkdown_2.14         cellranger_1.1.0       fastmatch_1.1-3
    ## [133] tidytree_0.3.9         rjson_0.2.21           nloptr_2.0.0
    ## [136] lifecycle_1.0.1        nlme_3.1-157           jsonlite_1.8.0
    ## [139] viridisLite_0.4.0      fansi_1.0.3            pillar_1.7.0
    ## [142] lattice_0.20-45        KEGGREST_1.36.0        fastmap_1.1.0
    ## [145] httr_1.4.2             GO.db_3.15.0           glue_1.6.2
    ## [148] png_0.1-7              iterators_1.0.14       bit_4.0.4
    ## [151] ggforce_0.3.3          stringi_1.7.6          sass_0.4.1
    ## [154] blob_1.2.3             memoise_2.0.1          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.