# signatureSearch: Environment for Gene Expression Searching Combined with Functional Enrichment Analysis

#### Authors: Yuzhu Duan, Dan S Evans, Nicholas J. Schork, Steven R. Cummings and Thomas Girke

#### Last update: 17 April, 2020

# 1 Introduction

## 1.1 Background

The signatureSearch package implements algorithms and data structures for performing gene expression signature (GES) searches, and subsequently interpreting the results functionally with specialized enrichment methods. These utilities are useful for studying the effects of genetic, chemical and environmental perturbations on biological systems. Specifically, in drug discovery they can be used for identifying novel modes of action (MOA) of bioactive compounds from reference databases such as LINCS containing the genome-wide GESs from tens of thousands of drug and genetic perturbations (Subramanian et al. 2017). A typical GES search (GESS) workflow can be divided into two major steps (Figure 1). First, GESS methods are used to identify perturbagens such as drugs that induce GESs similar to a query GES of interest. The queries can be drug-, disease- or phenotype-related GESs. Since the MOAs of most drugs in the corresponding reference databases are known, the resulting associations are useful to gain insights into pharmacological and/or disease mechanisms, and to develop novel drug repurposing approaches. Second, specialized functional enrichment analysis (FEA) methods using annotations systems, such as Gene Ontologies (GO), pathways or Disease Ontologies (DO), have been developed and implemented in this package to efficiently interpret GESS results. The latter are usually composed of lists of perturbagens (e.g. drugs) ranked by the similarity metric of the corresponding GESS method. Finally, network reconstruction functionalities are integrated for visualizing the final results, e.g. in form of drug-target networks. Figure 1 illustrates the major steps of a typical signature search workflow. For each GESS and FEA step, several alternative methods have been implemented in signatureSearch to allow users to choose the best possible workflow configuration for their research application. The individual search and enrichment methods are introduced in the corresponding sections of this vignette.

Figure 1: Overview of GESS and FEA methods. GES queries are used to search a perturbation-based GES reference database for perturbagens such as drugs inducing GESs similar to the query. To interpret the results, the GESS results are subjected to functional enrichment analysis (FEA) including drug set and target set enrichment analyses (DSEA, TSEA). Both identify functional categories (e.g. GO terms or KEGG pathways) over-represented in the GESS results. Subsequently, drug-target networks are reconstructed for visualization and interpretation.

## 1.2 Motivation and Design

Integrating the above described GESS and FEA methods into an R/Bioconductor package has several advantages. First, it simplifies the development of automated end-to-end workflows for conducting signature searches for many application areas. Second, it consolidates an extendable number of GESS and FEA algorithms into a single environment that allows users to choose the best selection of methods and parameter settings for a given research question. Third, the usage of generic data objects and classes improves maintainability and reproducibility of the provided functionalities, while the integration with the existing R/Bioconductor ecosystem maximizes their extensibility and reusability for other data analysis applications. Fourth, it provides access to several community perturbation reference databases along with options to build custom databases with support for most common gene expression profiling technologies (e.g. microarrays and RNA-Seq).

Figure 2 illustrates the design of the package with respect to its data containers and methods used by the individual signature search workflow steps. Briefly, expression profiles from genome-wide gene expression profiling technologies (e.g. RNA-Seq or microarrays) are used to build a reference database stored as HDF5 file. Commonly, a pre-built database can be used here that is provided by the associated signatureSearchData package. A search with a query signature against a reference database is initialized by declaring all parameter settings in a qSig search object. Users can choose here one of five different search algorithms implemented by signatureSearch. The nature of the query signature along with a chosen search method defines the type of expression data used for searching. For instance, a search with a query up/down gene set combined with the LINCS search method will be performed on sorted gene expression scores (see left panel in Figure 1). To minimize memory requirements and improve time performance, large reference databases are searched in batches with user-definable chunk sizes. The search results are stored in a gessResult object. The latter contains all information required to be processed by the downstream functional enrichment analysis (FEA) methods, here drug set and target set enrichment analysis (TSEA and DSEA) methods. The obtained functional enrichment results are stored as feaResult object that can be passed on to various drug-target network construction and visualization methods implemented in signatureSearch.

Figure 2: Design of signatureSearch. Gene expression profiles are stored in a reference database (here HDF5 file). Prebuilt databases for community GES collections, such as CMAP2 and LINCS, are provided by the affiliated signatureSearchData package. The GESS query parameters are defined in a qSig search object where users can choose among over five GESS methods (CMAP, LINCS, gCMAP, Fisher, and various correlation methods). Signature search results are stored in a gessResult object that can be functionally annotated with different TSEA (dup_hyperG, mGSEA, mabs) and DSEA (hyperG, GSEA) methods. The enrichment results are stored as `feaResult` object that can be used for drug target-networks analysis and visualization.

## 1.3 History of GES Databases

Lamb et al. (2006) generated one of the first GES databases called CMAP. Initially, it included GESs for 164 drugs screened against four mammalian cell lines (Lamb et al. 2006). A few years later CMAP was extended to CMAP2 containing GESs for 1,309 drugs and eight cell lines. More recently, a much larger GES database was released by the Library of Network-Based Cellular Signatures (LINCS) Consortium (Subramanian et al. 2017). In its initial release, the LINCS database contained perturbation-based GESs for 19,811 drugs tested on up to 70 cancer and non-cancer cell lines along with genetic perturbation experiments for several thousand genes. The number of compound dosages and time points considered in the assays has also been increased by 10-20 fold. The CMAP/CMAP2 databases use Affymetrix Gene Chips as expression platform. To scale from a few thousand to many hundred thousand GESs, the LINCS Consortium uses now the more economic L1000 assay. This bead-based technology is a low cost, high-throughput reduced representation expression profiling assay. It measures the expression of 978 landmark genes and 80 control genes by detecting fluorescent intensity of beads after capturing the ligation-mediated amplification products of mRNAs (Peck et al. 2006). The expression of 11,350 additional genes is imputed from the landmark genes by using as training data a collection of 12,063 Affymetrix gene chips (Edgar, Domrachev, and Lash 2002). The substantial scale-up of the LINCS project provides now many new opportunities to explore MOAs for a large number of known drugs and experimental drug-like small molecules.

## 1.4 Terminology

In the following text the term Gene Expression Signatures (GESs) can refer to at least four different situations of pre-processed gene expression data: (1) normalized gene expression intensity values (or counts for RNA-Seq); (2) \(log_2\) fold changes (LFC), z-scores or p-values obtained from analysis routines of differentially expressed genes (DEGs); (3) rank transformed versions of the expression values obtained under (1) and (2); and (4) gene identifier sets extracted from the top and lowest ranks under (3), such as `n` top up/down regulated DEGs. When unclear the text in this vignette specifies which of the four pre-processing states are used and/or the term GES refers to.

# 2 Getting Started

## 2.1 Package Install

As Bioconductor package `signatureSearch` can be installed with the `BiocManager::install()` function.

    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("signatureSearch")
    BiocManager::install("yduan004/signatureSearch", build_vignettes=TRUE)  # Installs from github

Next the package needs to be loaded into a user’s R session.

    library(signatureSearch)
    library(help = "signatureSearch")  # Lists package info
    vignette("signatureSearch")  # Opens vignette

    ## Warning: vignette 'signatureSearch' not found

## 2.2 Reference Databases

The helper package signatureSearchData provides access to pre-built GES databases, including CMAP2 and LINCS, that are stored on Bioconductor’s `ExperimentHub` as HDF5 files. Users can download these databases as follows.

    library(ExperimentHub); library(rhdf5)
    eh <- ExperimentHub()
    cmap <- eh[["EH3223"]]; cmap_expr <- eh[["EH3224"]]
    lincs <- eh[["EH3226"]]; lincs_expr <- eh[["EH3227"]]
    h5ls(lincs)

This will store the paths to the downloaded database files in the corresponding variables. The reference databases store the following information: (1) `lincs` contains moderated z-scores from differential expression (DE) analysis of 12,328 genes from 8,140 compound treatments of 30 cell lines corresponding to a total of 45,956 signatures; (2) `lincs_expr` contains gene expression intensity values from 5,925 compound treatments of 30 cell lines corresponding to a total of 38,824 signatures; (3) `cmap` contains \(log_2\) fold changes of 12,437 genes from 1,281 compound treatments of 5 cell lines corresponding to a total of 3,478 signatures; and (4) `cmap_expr` contains mean expression values from 1,309 drug treatments of 4 cell lines corresponding to a total of 3,587 signatures. To minimize redundancy in the `lincs` and `lincs_expr` databases, they were assembled from GESs corresponding to a compound dosage and treatment time of 10\(\mu\)M and 24h, respectively. If necessary one can create here easily database instances for all LINCS measurements. However, this will make the search results overwhelmingly complex which we wanted to avoid here. For details how the CMAP2 and LINCS databases were constructed, please refer to the vignette of the `signatureSearchData` package. The command `browseVignettes("signatureSearchData")` will open this vignette from R.

Custom databases can be built with the `build_custom_db` function. Here the user provides custom genome-wide gene expression data (e.g. for drug, disease or genetic perturbations) in a `data.frame` or `matrix`. The gene expression data can be most types of the pre-processed gene expression values described under section 1.4.

# 3 Signature Searches (GESS)

This tutorial section introduces the GESS algorithms implemented in signatureSearch. Currently, this includes five search algorithms, while additional ones will be added in the future. Based on the data types represented in the query and database, they can be classified into set- and correlation-based methods (see Figure 1). The first 4 methods described below are set-based, whereas the last one is a correlation-based method. We refer to a search method as set-based if at least one of the two data components (query and/or database) is composed of an identifier set (e.g. gene labels) that may be ranked or unranked. In contrast to this, correlation-based methods require quantitative values, usually of the same type such as normalized intensities, for both the query and the database entries. An advantage of the set-based methods is that their queries can be the highest and lowest ranking gene sets derived from a genome-wide profiling technology that may differ from the one used to generate the reference database. However, the precision of correlation methods often outperforms set-based methods. On the other hand, due to the nature of the expected input, correlation-based methods are usually only an option when both the query and database entries are based on the same or at least comparable technologies. In other words, set-based methods are more technology agnostic than correlation-based methods, but may not provide the best precision performance.

## 3.1 Test Query and Database

To minimize the run time of the following test code, a small toy database has been assembled from the LINCS database containing a total of 100 GESs from human SKB (muscle) cells. Of the 100 GESs in this toy database, 95 were random sampled and 5 were cherry-picked. The latter five are GESs from HDAC inhibitor treatments including the known drugs: vorinostat, rhamnetin, trichostatin A, pyroxamide, and HC toxin. To further reduce the size of the toy database, the number of its genes was reduced from 12,328 to 5,000 by random sampling. The query signature used in the sample code below is the vorinostat GES drawn from the toy database. For simplicity and to minimize the build time of this vignette, the following sample code uses a pre-generated instance of this toy database stored under the `extdata` directory of this package. The detailed code for generating the toy database and other custom instances of the LINCS database is given in the Supplementary Material section of this vignette.

The following imports the toy GES database into a `SummarizedExperiment` container (here `sample_db`). In addition, the test query set of up/down DEGs (here `upset` and `downset`) is extracted from the vorinostat GES entry in the database.

    db_path <- system.file("extdata", "sample_db.h5", package = "signatureSearch")
    # Load sample_db as `SummarizedExperiment` object
    library(SummarizedExperiment); library(HDF5Array)

    ## Loading required package: rhdf5

    sample_db <- SummarizedExperiment(HDF5Array(db_path, name="assay"))
    rownames(sample_db) <- HDF5Array(db_path, name="rownames")
    colnames(sample_db) <- HDF5Array(db_path, name="colnames")
    # get "vorinostat__SKB__trt_cp" signature drawn from toy database
    query_mat <- as.matrix(assay(sample_db[,"vorinostat__SKB__trt_cp"]))
    query = as.numeric(query_mat); names(query) = rownames(query_mat)
    upset <- head(names(query[order(-query)]), 150)
    downset <- tail(names(query[order(-query)]), 150)

## 3.2 CMAP Search Method

Lamb et al. (2006) introduced the gene expression-based search method known as Connectivity Map (CMap) where a GES database is searched with a query GES for similar entries (Lamb et al. 2006). Specifically, the GESS method from Lamb et al. (2006), here termed as CMAP, uses as query the two label sets of the most up- and down-regulated genes from a genome-wide expression experiment, while the reference database is composed of rank transformed expression profiles (e.g. ranks of LFC or z-scores). The actual GESS algorithm is based on a vectorized rank difference calculation. The resulting Connectivity Score expresses to what degree the query up/down gene sets are enriched on the top and bottom of the database entries, respectively. The search results are a list of perturbagens such as drugs that induce similar or opposing GESs as the query. Similar GESs suggest similar physiological effects of the corresponding perturbagens. As discussed in the introduction, these GES associations can be useful to uncover novel MOAs of drugs or treatments for diseases. Although several variants of the CMAP algorithm are available in other software packages including Bioconductor, the implementation provided by `signatureSearch` follows the original description of the authors as closely as possible. This allows to reproduce in our tests the search results from the corresponding CMAP2 web service of the Broad Institute.

In the following code block, the `qSig` function is used to generate a `qSig` object by defining (i) the query signature, (ii) the GESS method and (iii) the path to the reference database. Next, the query signature is used to search the reference database with the chosen GESS method. The type of the query signature and the reference database needs to meet the requirement of the search algorithm. In the chosen example the GESS method is CMAP, where the query signature needs to be the labels of up and down regulated genes, and the reference database contains rank transformed genome-wide expression profiles. Alternatively, the database can contain the genome-wide profiles themselves. In this case they will be transformed to gene ranks during the search. Note, in the given example the `db_path` variable stores the path to the toy database containing 100 GESs composed of z-scores.

    qsig_cmap <- qSig(query = list(upset=upset, downset=downset),
                      gess_method="CMAP", refdb=db_path)

    ## 150 / 150 genes in up set share identifiers with reference database

    ## 150 / 150 genes in down set share identifiers with reference database

    cmap <- gess_cmap(qSig=qsig_cmap, chunk_size=5000)
    cmap

    ## #
    ## # gessResult object
    ## #
    ## @result
    ## # A tibble: 100 x 9
    ##    pert   cell  type  trend raw_score scaled_score N_upset N_downset t_gn_sym
    ##    <chr>  <chr> <chr> <chr>     <dbl>        <dbl>   <int>     <int> <chr>
    ##  1 vorin~ SKB   trt_~ up        1.94         1         150       150 HDAC1; HDA~
    ##  2 resci~ SKB   trt_~ down     -0.295       -1         150       150 ACE
    ##  3 zuclo~ SKB   trt_~ down     -0.287       -0.972     150       150 ADRA1A; DR~
    ##  4 evoxi~ SKB   trt_~ down     -0.244       -0.828     150       150 <NA>
    ##  5 scoul~ SKB   trt_~ down     -0.242       -0.821     150       150 ADRA1D
    ##  6 gangl~ SKB   trt_~ down     -0.239       -0.809     150       150 <NA>
    ##  7 warfa~ SKB   trt_~ down     -0.234       -0.794     150       150 ARSE; CYP2~
    ##  8 trich~ SKB   trt_~ up        1.50         0.775     150       150 HDAC1; HDA~
    ##  9 endec~ SKB   trt_~ down     -0.224       -0.760     150       150 <NA>
    ## 10 HC-to~ SKB   trt_~ up        1.44         0.745     150       150 HDAC1
    ## # ... with 90 more rows
    ## @query    up gene set (150):      230 5357 2015 2542 1759 6195 80347 5050 5883 6251 ...
    ##           down gene set (150):    22864 9338 54793 10384 27000 10161 23363 114034 2997 23108 ...
    ##
    ## @gess_method      CMAP
    ##
    ## @refdb    C:/Users/biocbuild/bbs-3.11-bioc/tmpdir/Rtmp2brrwy/Rinst1f502fba781a/signatureSearch/extdata/sample_db.h5

    result(cmap)

    ## # A tibble: 100 x 9
    ##    pert   cell  type  trend raw_score scaled_score N_upset N_downset t_gn_sym
    ##    <chr>  <chr> <chr> <chr>     <dbl>        <dbl>   <int>     <int> <chr>
    ##  1 vorin~ SKB   trt_~ up        1.94         1         150       150 HDAC1; HDA~
    ##  2 resci~ SKB   trt_~ down     -0.295       -1         150       150 ACE
    ##  3 zuclo~ SKB   trt_~ down     -0.287       -0.972     150       150 ADRA1A; DR~
    ##  4 evoxi~ SKB   trt_~ down     -0.244       -0.828     150       150 <NA>
    ##  5 scoul~ SKB   trt_~ down     -0.242       -0.821     150       150 ADRA1D
    ##  6 gangl~ SKB   trt_~ down     -0.239       -0.809     150       150 <NA>
    ##  7 warfa~ SKB   trt_~ down     -0.234       -0.794     150       150 ARSE; CYP2~
    ##  8 trich~ SKB   trt_~ up        1.50         0.775     150       150 HDAC1; HDA~
    ##  9 endec~ SKB   trt_~ down     -0.224       -0.760     150       150 <NA>
    ## 10 HC-to~ SKB   trt_~ up        1.44         0.745     150       150 HDAC1
    ## # ... with 90 more rows

The search result is stored in a `gessResult` object (here named `cmap`) containing the following components: search result table, query signature, name of the GESS method and path to the reference database. The `result` accessor function can be used to extract the search result table from the `gessResult` object. This table contains the search results for each perturbagen (here drugs) in the reference database ranked by their signature similarity to the query. For the CMAP method, the similarity metrics are `raw_score` and `scaled_score`. The raw score represents the bi-directional enrichment score (Kolmogorov-Smirnov statistic) for a given up/down query signature. Under the `scaled_score` column, the `raw_score` has been scaled to values from 1 to -1 by dividing positive scores and negative scores with the maximum positive score and the absolute value of the minimum negative score, respectively. The remaining columns in the search result table contain the following information. `pert`: name of perturbagen (e.g. drug) in the reference database; `cell`: acronym of cell type; `type`: perturbation type, e.g. compound treatment is `trt_cp`; `trend`: up or down when reference signature is positively or negatively connected with the query signature, respectively; `N_upset` or `N_downset`: number of genes in the query up or down sets, respectively; `t_gn_sym`: gene symbols of the corresponding drug targets.

## 3.3 LINCS Search Method

Subramanian et al. (2017) introduced a more complex GESS algorithm, here referred to as LINCS. While related to CMAP, there are several important differences among the two approaches. First, LINCS weights the query genes based on the corresponding differential expression scores of the GESs in the reference database (e.g. LFC or z-scores). Thus, the reference database used by LINCS needs to store the actual score values rather than their ranks. Another relevant difference is that the LINCS algorithm uses a bi-directional weighted Kolmogorov-Smirnov enrichment statistic (ES) as similarity metric. To the best of our knowledge, the `LINCS` search functionality in `signatureSearch` provides the first downloadable standalone software implementation of this algorithm.

In the following example the `qSig` object for the LINCS method is initialized the same way as the corresponding CMAP object above with the exception that “LINCS” needs to be specified under the `gess_method` argument.

    qsig_lincs <- qSig(query = list(upset=upset, downset=downset),
                       gess_method="LINCS", refdb=db_path)

    ## 150 / 150 genes in up set share identifiers with reference database

    ## 150 / 150 genes in down set share identifiers with reference database

    lincs <- gess_lincs(qsig_lincs, sortby="NCS", tau=FALSE)
    result(lincs)

    ## # A tibble: 100 x 12
    ##    pert  cell  type  trend   WTCS WTCS_Pval WTCS_FDR   NCS NCSct N_upset
    ##    <chr> <chr> <chr> <chr>  <dbl>     <dbl>    <dbl> <dbl> <dbl>   <int>
    ##  1 vori~ SKB   trt_~ up     1       0.      0         2.57  2.57     150
    ##  2 tric~ SKB   trt_~ up     0.863   5.96e-6 0.000180  2.22  2.22     150
    ##  3 HC-t~ SKB   trt_~ up     0.857   6.22e-6 0.000180  2.20  2.20     150
    ##  4 pyro~ SKB   trt_~ up     0.625   1.63e-5 0.000180  1.60  1.60     150
    ##  5 zucl~ SKB   trt_~ down  -0.321   8.29e-5 0.000319 -1.19 -1.19     150
    ##  6 resc~ SKB   trt_~ down  -0.319   9.60e-5 0.000356 -1.18 -1.18     150
    ##  7 APHA~ SKB   trt_~ up     0.445   2.42e-5 0.000180  1.14  1.14     150
    ##  8 epot~ SKB   trt_~ down  -0.308   2.40e-4 0.000801 -1.14 -1.14     150
    ##  9 scop~ SKB   trt_~ up     0.416   2.54e-5 0.000180  1.07  1.07     150
    ## 10 I-07~ SKB   trt_~ up     0.408   2.58e-5 0.000180  1.05  1.05     150
    ## # ... with 90 more rows, and 2 more variables: N_downset <int>, t_gn_sym <chr>

The search results are stored in a `gessResult` object as under the CMAP example above. The similarity scores stored in the LINCS result table are summarized here. `WTCS`: Weighted Connectivity Score; `WTCS_Pval`: nominal p-value of WTCS; `WTCS_FDR`: false discovery rate of `WTCS_Pval`; `NCS`: normalized connectivity score; `NCSct`: NCS summarized across cell types; `Tau`: enrichment score standardized for a given database. The latter is only included in the result table if `tau=TRUE` in a `gess_lincs` function call. The example given is run with `tau=FALSE`, because the tau values are only meaningful when the complete LINCS database is used which is not the case for the toy database.

The following provides a more detailed description of the similarity scores computed by the LINCS method. Additional details are available in the Supplementary Material Section of the Subramanian et al. (2017) paper.

`WTCS`: The Weighted Connectivity Score is a bi-directional ES for an up/down query set. If the ES values of an up set and a down set are of different signs, then WTCS is (ESup-ESdown)/2, otherwise, it is 0. WTCS values range from -1 to 1. They are positive or negative for signatures that are positively or inversely related, respectively, and close to zero for signatures that are unrelated.

`WTCS_Pval` and `WTCS_FDR`: The nominal p-value of the WTCS and the corresponding false discovery rate (FDR) are computed by comparing the WTCS against a null distribution of WTCS values obtained from a large number of random queries (e.g. 1000).

`NCS`: To make connectivity scores comparable across cell types and perturbation types, the scores are normalized. Given a vector of \(WTCS\) values \(w\) resulting from a query, the values are normalized within each cell line \(c\) and perturbagen type \(t\) to obtain the Normalized Connectivity Score (\(NCS\)) by dividing the \(WTCS\) value by the signed mean of the \(WTCS\) values within the subset of signatures in the reference database corresponding to \(c\) and \(t\).

`NCSct`: The NCS is summarized across cell types as follows. Given a vector of \(NCS\) values for perturbagen \(p\), relative to query \(q\), across all cell lines \(c\) in which \(p\) was profiled, a cell-summarized connectivity score is obtained using a maximum quantile statistic. It compares the 67 and 33 quantiles of \(NCSp,c\) and retains whichever is of higher absolute magnitude.

`Tau`: The standardized score Tau compares an observed \(NCS\) to a large set of \(NCS\) values that have been pre-computed for a specific reference database. The query results are scored with Tau as a standardized measure ranging from 100 to -100. A Tau of 90 indicates that only 10% of reference perturbations exhibit stronger connectivity to the query. This way one can make more meaningful comparisons across query results.

## 3.4 gCMAP Search Method

The Bioconductor gCMAP (Sandmann et al. 2014) package provides access to a related but not identical implementation of the original CMAP algorithm proposed by Lamb et al. (2006). It uses as query a rank transformed GES and the reference database is composed of the labels of up and down regulated DEG sets. This is the opposite situation of the CMAP method, where the query is composed of the labels of up and down regulated DEGs and the database contains rank transformed GESs.

In case of the gCMAP GESS method, the query GES is a matrix with a single column representing gene ranks from a biological state of interest, here vorinostat treatment in SKB cells. The corresponding gene labels are stored in the row name slot of the matrix. Instead of ranks one can provide scores (e.g. z-scores) as in the example given below. In such a case the scores will be internally transformed to ranks. The reference database consists of gene label sets that were extracted from the toy databases by applying a `higher` and `lower` filter, here set to `1` and `-1`, respectively.

    qsig_gcmap <- qSig(query = query_mat, gess_method = "gCMAP", refdb = db_path)
    gcmap <- gess_gcmap(qsig_gcmap, higher=1, lower=-1)
    result(gcmap)

    ## # A tibble: 100 x 9
    ##    pert     cell  type  trend effect  nSet nFound signed t_gn_sym
    ##    <chr>    <chr> <chr> <chr>  <dbl> <dbl>  <dbl> <lgl>  <chr>
    ##  1 vorinos~ SKB   trt_~ up     1      1073   1073 TRUE   HDAC1; HDAC10; HDAC11;~
    ##  2 warfarin SKB   trt_~ down  -1       105    105 TRUE   ARSE; CYP2C8; CYP2C9; ~
    ##  3 D-609    SKB   trt_~ down  -0.720   121    121 TRUE   <NA>
    ##  4 TC-2559  SKB   trt_~ down  -0.681   367    367 TRUE   CHRNA4
    ##  5 rescinn~ SKB   trt_~ down  -0.618   707    707 TRUE   ACE
    ##  6 trichos~ SKB   trt_~ up     0.578  1408   1408 TRUE   HDAC1; HDAC10; HDAC2; ~
    ##  7 amiodar~ SKB   trt_~ down  -0.562   686    686 TRUE   ABCG2; ADRB1; CACNA1H;~
    ##  8 zuclope~ SKB   trt_~ down  -0.545   464    464 TRUE   ADRA1A; DRD1; HRH1
    ##  9 progest~ SKB   trt_~ down  -0.535   235    235 TRUE   AKR1C1; AKR1C2; AKR1D1~
    ## 10 ganglio~ SKB   trt_~ down  -0.523   587    587 TRUE   <NA>
    ## # ... with 90 more rows

As in the other search methods, the gCMAP results are stored in a `gessResult` object. The columns in the corresponding search result table, that are specific to the gCMAP method, contain the following information. `effect`: scaled bi-directional enrichment score corresponding to the `scaled_score` under the CMAP result; `nSet`: number of genes in the reference gene sets after applying the higher and lower cutoff; `nFound`: number of genes in the reference gene sets that are present in the query signature; `signed`: whether the gene sets in the reference database have signs, e.g. representing up and down regulated genes when computing scores.

## 3.5 Fisher Search Method

Fisher’s exact test (Graham J. G. Upton 1992) can also be used to search a GES databases for entries that are similar to a query GES. In this case both the query and the database are composed of gene label sets, such as DEG sets.

In the following example, both the `query` and the `refdb` used under the `qSig` call are genome-wide signatures, here z-scores. The actual gene sets required for the Fisher’s exact test are obtained by setting the `higher` and `lower` cutoffs to 1 and -1, respectively.

    qsig_fisher <- qSig(query = query_mat, gess_method = "Fisher", refdb = db_path)
    fisher <- gess_fisher(qSig=qsig_fisher, higher=1, lower=-1)
    result(fisher)

    ## # A tibble: 100 x 12
    ##    pert  cell  type  trend      pval      padj effect     LOR  nSet nFound
    ##    <chr> <chr> <chr> <chr>     <dbl>     <dbl>  <dbl>   <dbl> <dbl>  <dbl>
    ##  1 vori~ SKB   trt_~ over  0.        0.         37.5  Inf      1073   1073
    ##  2 tric~ SKB   trt_~ over  8.03e-181 4.01e-179  28.7    2.10   1408    693
    ##  3 HC-t~ SKB   trt_~ over  2.09e-136 6.96e-135  24.9    1.78   1786    735
    ##  4 SPB0~ SKB   trt_~ under 3.92e- 26 9.79e- 25 -10.6   -1.12    930     88
    ##  5 MDL-~ SKB   trt_~ under 9.40e- 24 1.88e- 22 -10.0   -0.770  1788    247
    ##  6 clof~ SKB   trt_~ under 2.81e- 18 4.68e- 17  -8.72  -0.648  1912    289
    ##  7 pyro~ SKB   trt_~ over  1.24e- 17 1.78e- 16   8.55   0.813   623    220
    ##  8 I-07~ SKB   trt_~ over  1.56e- 12 1.95e- 11   7.07   0.545  1207    349
    ##  9 TUL-~ SKB   trt_~ under 1.01e- 11 1.12e- 10  -6.81  -0.694   855    112
    ## 10 form~ SKB   trt_~ under 3.19e- 11 3.19e- 10  -6.64  -0.470  2210    379
    ## # ... with 90 more rows, and 2 more variables: signed <lgl>, t_gn_sym <chr>

The columns in the result table specific to the Fisher method include the following information. `pval`: p-value of the Fisher’s exact test; `padj`: p-value adjusted for multiple hypothesis testing using R’s `p.adjust` function with the Benjamini & Hochberg (BH) method; `effect`: z-score based on the standard normal distribution; `LOR`: log odds ratio.

If the `query` contains the labels of up and down regulated genes then the two sets can be provided as a list. Internally, they will be combined into a single unsigned set, while the reference database is processed the same way as in the previous example.

    qsig_fisher2 <- qSig(query = list(upset=upset, downset=downset),
                         gess_method = "Fisher", refdb = db_path)
    fisher2 <- gess_fisher(qSig=qsig_fisher2, higher=1, lower=-1)
    result(fisher2)

## 3.6 Correlation-based Search Method

Correlation-based similarity metrics, such as Spearman or Pearson coefficients, can also be used as GESS methods. As non-set-based methods, they require quantitative gene expression values for both the query and the database entries, that usually need to be of the same type to obtain meaningful results, such as normalized intensities or read counts from microarrays or RNA-Seq experiments, respectively. For correlation searches to work, it is important that both the query and reference database contain the same type of gene identifiers. The expected data structure of the query is a matrix with a single numeric column and the gene labels (e.g. Entrez Gene IDs) in the row name slot. For convenience, the correlation-based searches can either be performed with the full set of genes represented in the database or a subset of them. The latter can be useful to focus the computation for the correlation values on certain genes of interest such as a DEG set or the genes in a pathway of interest. For comparing the performance of different GESS methods, it can also be advantageous to subset the genes used for a correlation-based search to same set used in a set-based search, such as the up/down DEGs used in a LINCS GESS. This way the search results of correlation- and set-based methods can be more comparable because both are provided with equivalent information content.

### 3.6.1 CORall

The following example runs a correlation-based search with the Spearman method using all genes present in the reference database. The GESs used for both the `query` and `refdb` are z-scores. The `gess_cor` function also supports Pearson and Kendall correlation coefficients by assigning the corresponding names to the `method` argument. For details, users want to consult the help file of the `gess_cor` function.

    qsig_sp <- qSig(query = query_mat, gess_method = "Cor", refdb = db_path)
    sp <- gess_cor(qSig=qsig_sp, method="spearman")
    result(sp)

    ## # A tibble: 100 x 6
    ##    pert        cell  type   trend cor_score t_gn_sym
    ##    <chr>       <chr> <chr>  <chr>     <dbl> <chr>
    ##  1 vorinostat  SKB   trt_cp up        1     HDAC1; HDAC10; HDAC11; HDAC2; HDAC3~
    ##  2 trichostat~ SKB   trt_cp up        0.702 HDAC1; HDAC10; HDAC2; HDAC3; HDAC4;~
    ##  3 HC-toxin    SKB   trt_cp up        0.634 HDAC1
    ##  4 pyroxamide  SKB   trt_cp up        0.325 HDAC1; HDAC3; HDAC6
    ##  5 APHA-compo~ SKB   trt_cp up        0.152 HDAC8
    ##  6 benzonatate SKB   trt_cp up        0.151 SCN5A
    ##  7 rescinnami~ SKB   trt_cp down     -0.136 ACE
    ##  8 evoxine     SKB   trt_cp down     -0.134 <NA>
    ##  9 scopolamin~ SKB   trt_cp up        0.125 <NA>
    ## 10 warfarin    SKB   trt_cp down     -0.123 ARSE; CYP2C8; CYP2C9; CYP4F2; F3; F~
    ## # ... with 90 more rows

The column specific to the correlation-based search methods contains the following information. `cor_score`: correlation coefficient based on the method defined in the `gess_cor` function call.

### 3.6.2 CORsub

To perform a correlation-based search on a subset of genes represented in the database, one can simply provide the chosen gene subset in the query. During the search the database entries will be subsetted to the genes provided in the query signature. The given example uses a query GES that is subsetted to the genes with the 150 highest and lowest z-scores.

    # Subset z-scores of 150 up and down gene sets from
    # "vorinostat__SKB__trt_cp" signature.
    query_mat_sub <- as.matrix(query_mat[c(upset, downset),])
    qsig_spsub <- qSig(query = query_mat_sub, gess_method = "Cor", refdb = db_path)
    spsub <- gess_cor(qSig=qsig_spsub, method="spearman")
    result(spsub)

    ## # A tibble: 100 x 6
    ##    pert        cell  type   trend cor_score t_gn_sym
    ##    <chr>       <chr> <chr>  <chr>     <dbl> <chr>
    ##  1 vorinostat  SKB   trt_cp up        1     HDAC1; HDAC10; HDAC11; HDAC2; HDAC3~
    ##  2 trichostat~ SKB   trt_cp up        0.885 HDAC1; HDAC10; HDAC2; HDAC3; HDAC4;~
    ##  3 HC-toxin    SKB   trt_cp up        0.856 HDAC1
    ##  4 pyroxamide  SKB   trt_cp up        0.646 HDAC1; HDAC3; HDAC6
    ##  5 APHA-compo~ SKB   trt_cp up        0.442 HDAC8
    ##  6 benzonatate SKB   trt_cp up        0.333 SCN5A
    ##  7 MD-040      SKB   trt_cp up        0.307 <NA>
    ##  8 K784-3187   SKB   trt_cp up        0.297 <NA>
    ##  9 scopolamin~ SKB   trt_cp up        0.297 <NA>
    ## 10 fasudil     SKB   trt_cp up        0.288 CDC42BPB; MYLK3; PKIA; PRKACA; PRKG~
    ## # ... with 90 more rows

## 3.7 Summary of Search Results

Although the toy database is artificially small, one can use the above search results for a preliminary performance assessment of the different GESS methods in ranking drugs based on known modes of action (MOA). Four of the five cherry-picked HDAC inhibitors (vorinostat, trichostatin-a, HC-toxin, pyroxamide) were ranked among the top 10 ranking drugs in the search results of the LINCS, Fisher and Spearman correlation methods. If generalizable, this result implies a promising performance of these search methods for grouping drugs by their MOA categories. In addition, the LINCS and Spearman methods were able to rank another HDAC inhibitor, APHA-compound-8, at the top of their search results, indicating a better sensitivity of these two methods compared to the other methods.

## 3.8 GESS Result Visualization

The `gess_res_vis` function allows to summarize the ranking scores of selected perturbagens for GESS results across cell types along with cell type classifications, such as normal and tumor cells. In the following plot (Figure 3) the perturbagens are drugs (along x-axis) and the ranking scores are LINCS’ NCS values (y-axis). For each drug the NCS values are plotted for each cell type as differently colored dots, while their shape indicates the cell type class. Note, the code for generating the plot is not evaluated here since the toy database used by this vignette contains only treatments for one cell type (here SKB cells). This would result in a not very informative plot. To illustrate the full potential of the `gess_res_vis` function, the following code section applies to a search where the `vorinostat` signature was used to query with the `gess_lincs` method the full LINCS database. Subsequently, the search result is processed by the `gess_res_vis` function to generate the plot shown in Figure 3.

    vor_qsig_full <- qSig(query = list(upset=upset, downset=downset),
                       gess_method="LINCS", refdb="lincs")
    vori_res_full <- gess_lincs(qSig=vor_qsig_full, sortby="NCS", tau=TRUE)
    vori_tb <- result(vori_res_full)
    drugs_top10 <- unique(result(lincs)$pert)[1:10]
    drugs_hdac <- c("panobinostat","mocetinostat","ISOX","scriptaid","entinostat",
          "belinostat","HDAC3-selective","tubastatin-a","tacedinaline","depudecin")
    drugs = c(drugs_top10, drugs_hdac)
    gess_res_vis(vori_tb, drugs = drugs, col = "NCS")

Figure 3: Summary of NCS scores across cell types (y-axis) for selected drugs (x-axis). The plot shown here is based on a search result where the vorinostat signature was used to query the entire LINCS database with the `gess_lincs` method. The drugs included here are the 10 top ranking drugs of the search result plus 10 cherry-picked drugs that are all HDAC inhibitors. Additional details are provided in the text of this sub-section.

## 3.9 Batch Processing of GESSs

Multiple GES queries can be processed sequentially or in parallel mode. Parallel evaluations can substantially reduce processing times. The parallelization techniques covered in this vignette, are based on utilities of the `BiocParallel` and `batchtools` packages. For demonstration purposes the following example uses a small batch query containing several GESs. First, this batch query is processed sequentially without any parallelization using a simple `lapply` loop. Next, the same query is processed in parallel mode using multiple CPU cores of a single machine. The third option demonstrates how this query can be processed in parallel mode on multiple machines of a computer cluster with a workload management (queueing) system (e.g. such as Slurm or Torque).

### 3.9.1 Sequential Processing

The following processes a small toy batch query with the `LINCS` method sequentially in an `lapply` loop. The object `batch_queries` is a list containing the two sample GESs `q1` and `q2` that are composed of entrez gene identifiers. The search results are written to tab-delimited tabular files under a directory called `batch_res`. The name and path of this directory can be changed as needed.

    library(readr)
    batch_queries <- list(q1=list(upset=c("23645", "5290"), downset=c("54957", "2767")),
                          q2=list(upset=c("27101","65083"), downset=c("5813", "84")))
    refdb <- system.file("extdata", "sample_db.h5", package="signatureSearch")
    gess_list <- lapply(seq_along(batch_queries), function(i){
        qsig_lincs <- qSig(query = batch_queries[[i]],
                       gess_method="LINCS", refdb=refdb)
        lincs <- gess_lincs(qsig_lincs, sortby="NCS", tau=TRUE)
        if(!dir.exists("batch_res")){
            dir.create("batch_res")
        }
        write_tsv(result(lincs), paste0("batch_res/lincs_res_", i, ".tsv"))
        return(result(lincs))
    })

### 3.9.2 Parallelization with Multiple CPU Cores

The GESSs from the previous example can be accelerated by taking advantage of multiple CPU cores available on a single computer system. The parallel evaluation happens in the below `bplapply` loop defined by the `BiocParallel` package. For this approach, all processing instructions are encapsulated in a function named `f_bp` that will be executed in the `bplapply` loop. As before, the search results are written to tab-delimited tabular files under a directory called `batch_res`. The name and path of this directory can be changed as needed. For more background information on this and the following parallelization options, users want to consult the vignette of the `BiocParallel` package.

    library(BiocParallel)
    f_bp <- function(i){
        qsig_lincs <- qSig(query = batch_queries[[i]],
                       gess_method="LINCS", refdb=refdb)
        lincs <- gess_lincs(qsig_lincs, sortby="NCS", tau=TRUE)
        if(!dir.exists("batch_res")){
            dir.create("batch_res")
        }
        write_tsv(result(lincs), paste0("batch_res/lincs_res_", i, ".tsv"))
        return(result(lincs))
    }
    gess_list <- bplapply(seq_along(batch_queries), f_bp, BPPARAM = MulticoreParam(workers = 2))

### 3.9.3 Parallelization with Multiple Computer Systems

In addition to utilizing multiple CPU cores of a single machine, one can further accelerate the processing by taking advantage of multiple computer systems (nodes) available on a computer cluster, where a queueing systems takes care of the load balancing.

In the following example, `Njobs` sets the number of independent processes to be run on the cluster, and `ncpus` defines the number of CPU cores to be used by each process. The chosen example will run 2 processes each utilizing 2 CPU cores. If `batch_queries` contains sufficient GESs and the corresponding computing resources are available on a cluster, then the given example process will utilize in total 4 CPU cores. Note, the given sample code will work on most queueing systems as it is based on utilities from the `batchtools` package. The latter supports template files (`*.tmpl`) for defining the run parameters of different schedulers. To run this code, one needs to have both a `conf` file (see `.batchtools.conf.R` samples here) and a `template` file (see `*.tmpl` samples here) for the queueing system available on a cluster. The following example uses the sample `conf` and `template` files for the Slurm scheduler.

For additional details on parallelizing computations on clusters, users want to consult the vignettes of the `batchtools` and `BiocParallel` packages.

    library(batchtools)
    batch_queries <- list(q1=list(upset=c("23645", "5290"), downset=c("54957", "2767")),
                          q2=list(upset=c("27101","65083"), downset=c("5813", "84")))
    refdb <- system.file("extdata", "sample_db.h5", package="signatureSearch")

    f_bt <- function(i, batch_queries, refdb){
        library(signatureSearch)
        library(readr)
        qsig_lincs <- qSig(query = batch_queries[[i]],
                       gess_method="LINCS", refdb=refdb)
        lincs <- gess_lincs(qsig_lincs, sortby="NCS", tau=TRUE)
        if(!dir.exists("batch_res")){
            dir.create("batch_res")
        }
        write_tsv(result(lincs), paste0("batch_res/lincs_res_", i, ".tsv"))
        return(result(lincs)) # or return()
    }

Copy the `conf` and `template` files for Slurm to current working directory.

    file.copy(system.file("extdata", ".batchtools.conf.R", package="signatureSearch"), ".")
    file.copy(system.file("extdata", "slurm.tmpl", package="signatureSearch"), ".")

Create the registry and submit jobs.

    reg <- makeRegistry(file.dir="reg_batch", conf.file=".batchtools.conf.R")
    # reg <- loadRegistry(file.dir="reg_batch", conf.file=".batchtools.conf.R", writeable=TRUE)
    Njobs <- 1:2
    ids <- batchMap(fun=f_bt, Njobs, more.args = list(
              batch_queries=batch_queries, refdb=refdb))
    done <- submitJobs(ids, reg=reg, resources=list(
                partition="intel", walltime=120, ntasks=1, ncpus=2, memory=10240))
    getStatus()
    waitForJobs() # Wait until all jobs are completed
    res1 <- loadResult(1)
    unlink(c(".batchtools.conf.R", "slurm.tmpl"))

# 4 Functional Enrichment (FEA)

GESS results are lists of perturbagens (here drugs) ranked by their signature similarity to a query signature of interest. Interpreting these search results with respect to the cellular networks and pathways affected by the top ranking drugs is difficult. To overcome this challenge, the knowledge of the target proteins of the top ranking drugs can be used to perform functional enrichment analysis (FEA) based on community annotation systems, such as Gene Ontologies (GO), pathways (e.g. KEGG, Reactome), drug MOAs or Pfam domains. For this, the ranked drug sets are converted into target gene/protein sets to perform Target Set Enrichment Analysis (TSEA) based on a chosen annotation system. Alternatively, the functional annotation categories of the targets can be assigned to the drugs directly to perform Drug Set Enrichment Analysis (DSEA). Although TSEA and DSEA are related, their enrichment results can be distinct. This is mainly due to duplicated targets present in the test sets of the TSEA methods, whereas the drugs in the test sets of DSEA are usually unique. Additional reasons include differences in the universe sizes used for TSEA and DSEA.

Importantly, the duplications in the test sets of the TSEA are due to the fact that many drugs share the same target proteins. Standard enrichment methods would eliminate these duplications since they assume uniqueness in the test sets. Removing duplications in TSEA would be inappropriate since it would erase one of the most important pieces of information of this approach. To solve this problem, we have developed and implemented in the `signatureSearch` package a weighting method for duplicated targets, where the weighting is proportional to the frequency of the targets in the test set.

To perform TSEA and DSEA, drug-target annotations are essential. They can be obtained from several sources, including DrugBank, ChEMBL, STITCH, and the Touchstone dataset from the LINCS project (Wishart et al. 2018; Gaulton et al. 2017; Kuhn et al. 2010; Subramanian et al. 2017). Most drug-target annotations provide UniProt identifiers for the target proteins. They can be mapped, if necessary via their encoding genes, to the chosen functional annotation categories, such as GO or KEGG. To minimize bias in TSEA or DSEA, often caused by promiscuous binders, it can be beneficial to remove drugs or targets that bind to large numbers of distinct proteins or drugs, respectively. To conduct TSEA and DSEA efficiently, `signatureSearch` and its helper package `signatureSearchData` provide several convenience utilities along with drug-target lookup resources for automating the mapping from drug sets to target sets to functional categories.

Note, most FEA tests involving proteins in their test sets are performed on the gene level in `signatureSearch`. This way one can avoid additional duplications due to many-to-one relationships among proteins and their encoding genes. For this, the corresponding functions in `signatureSearch` will usually translate target protein sets into their encoding gene sets using identifier mapping resources from R/Bioconductor, such as the `org.Hs.eg.db` annotation package. Because of this as well as simplicity, the text in the vignette and help files of this package will refer to the targets of drugs almost interchangeably as proteins or genes, even though the former are the direct targets and the latter only the indirect targets of drugs.

## 4.1 TSEA

The following introduces how to perform TSEA on drug-based GESS results using as functional annotation systems GO and KEGG pathways. The enrichment tests are performed with three widely used algorithms that have been modified in `signatureSearch` to take advantage of duplication information present in the test sets used for TSEA. The relevance of these target duplications is explained above. The specialized enrichment algorithms include Duplication Adjusted Hypergeometric Test (`dup_hyperG`), Modified Gene Set Enrichment Analysis (`mGSEA`) and MeanAbs (`mabs`).

### 4.1.1 Hypergeometric Test

The classical hypergeometric test assumes uniqueness in its gene/protein test sets. Its p-value is calculated according to equation

\[\begin{equation} p=\sum_{k=x}^{n} \frac {{D \choose k}{{N-D} \choose {n-k}}}{{N \choose n}}. \end{equation}\]

In case of GO term enrichment analysis the individual variables are assigned the following components. \(N\) is the total number of genes/proteins contained in the entire annotation universe; \(D\) is the number of genes annotated at a specific GO node; \(n\) is the total number of genes in the test set; and \(x\) is the number of genes in the test set annotated at a specific GO node. To maintain the duplication information in the test sets used for TSEA, the values of \(n\) and \(x\) in the above equation are adjusted by the frequency of the target proteins in the test set. Effectively, the approach removes the duplications, but maintains their frequency information in form of weighting values.

#### 4.1.1.1 With GO

The following example code uses the \(n\) top ranking drugs (here \(n=10\)) from the LINCS GESS result as input to the `tsea_dup_hyperG` method. Internally, the latter converts the drug set to a target set, and then computes for it enrichment scores for each MF GO term based on the hypergeometric distribution. The enrichment results are stored in a `feaResult` object. It contains the organism information of the annotation system, and the ontology type of the GO annotation system. If the annotation system is KEGG, the latter will be “KEGG”. The object also stores the input drugs used for the enrichment test, as well as their target information.

    drugs <- unique(result(lincs)$pert[1:10])
    dup_hyperG_res <- tsea_dup_hyperG(drugs = drugs, universe = "Default",
                                      type = "GO", ont="MF", pvalueCutoff=0.05,
                                      pAdjustMethod="BH", qvalueCutoff = 0.1,
                                      minGSSize = 10, maxGSSize = 500)

    ## Loading required package: org.Hs.eg.db

    ## Loading required package: AnnotationDbi

    ##

    dup_hyperG_res

    ## #
    ## # Functional Enrichment Analysis
    ## #
    ## #...@organism     Homo sapiens
    ## #...@ontology     MF
    ## #...@drugs    chr [1:10] "vorinostat" "trichostatin-a" "hc-toxin" "pyroxamide" ...
    ## #...@targets      chr [1:41] "HDAC8" "TUBA1A" "TUBA1B" "TUBA1C" "TUBA3C" "TUBA4A" "TUBA8" ...
    ## #...71 enriched terms found
    ## # A tibble: 71 x 10
    ##    ont   ID    Description GeneRatio BgRatio   pvalue p.adjust   qvalue itemID
    ##    <chr> <chr> <chr>       <chr>     <chr>      <dbl>    <dbl>    <dbl> <chr>
    ##  1 MF    GO:0~ NAD-depend~ 26/41     16/193~ 0.       0.       0.       HDAC8~
    ##  2 MF    GO:0~ histone de~ 26/41     11/193~ 0.       0.       0.       HDAC8~
    ##  3 MF    GO:0~ NAD-depend~ 26/41     11/193~ 0.       0.       0.       HDAC8~
    ##  4 MF    GO:0~ NAD-depend~ 26/41     17/193~ 0.       0.       0.       HDAC8~
    ##  5 MF    GO:0~ histone de~ 26/41     44/193~ 9.55e-63 1.65e-60 5.49e-61 HDAC8~
    ##  6 MF    GO:0~ protein de~ 26/41     45/193~ 2.26e-62 3.41e-60 1.14e-60 HDAC8~
    ##  7 MF    GO:0~ deacetylas~ 26/41     58/193~ 2.03e-58 2.72e-56 9.08e-57 HDAC8~
    ##  8 MF    GO:0~ hydrolase ~ 26/41     89/193~ 1.85e-52 1.86e-50 6.20e-51 HDAC8~
    ##  9 MF    GO:0~ hydrolase ~ 26/41     147/19~ 4.54e-46 3.42e-44 1.14e-44 HDAC8~
    ## 10 MF    GO:0~ histone de~ 20/41     112/19~ 7.39e-35 5.24e-33 1.75e-33 HDAC1~
    ## # ... with 61 more rows, and 1 more variable: Count <int>

The `result` accessor function can be used to extract a tabular result from the `feaResult` object. The rows of this result table contain the functional categories (e.g. GO terms or KEGG pathways) ranked by the corresponding enrichment statistic.

    result(dup_hyperG_res)

    ## # A tibble: 71 x 10
    ##    ont   ID    Description GeneRatio BgRatio   pvalue p.adjust   qvalue itemID
    ##    <chr> <chr> <chr>       <chr>     <chr>      <dbl>    <dbl>    <dbl> <chr>
    ##  1 MF    GO:0~ NAD-depend~ 26/41     16/193~ 0.       0.       0.       HDAC8~
    ##  2 MF    GO:0~ histone de~ 26/41     11/193~ 0.       0.       0.       HDAC8~
    ##  3 MF    GO:0~ NAD-depend~ 26/41     11/193~ 0.       0.       0.       HDAC8~
    ##  4 MF    GO:0~ NAD-depend~ 26/41     17/193~ 0.       0.       0.       HDAC8~
    ##  5 MF    GO:0~ histone de~ 26/41     44/193~ 9.55e-63 1.65e-60 5.49e-61 HDAC8~
    ##  6 MF    GO:0~ protein de~ 26/41     45/193~ 2.26e-62 3.41e-60 1.14e-60 HDAC8~
    ##  7 MF    GO:0~ deacetylas~ 26/41     58/193~ 2.03e-58 2.72e-56 9.08e-57 HDAC8~
    ##  8 MF    GO:0~ hydrolase ~ 26/41     89/193~ 1.85e-52 1.86e-50 6.20e-51 HDAC8~
    ##  9 MF    GO:0~ hydrolase ~ 26/41     147/19~ 4.54e-46 3.42e-44 1.14e-44 HDAC8~
    ## 10 MF    GO:0~ histone de~ 20/41     112/19~ 7.39e-35 5.24e-33 1.75e-33 HDAC1~
    ## # ... with 61 more rows, and 1 more variable: Count <int>

The columns in the result table, extracted from the `feaResult` object, contain the following information. Note, some columns are only present in the result tables of specific FEA methods. `ont`: in case of GO one of BP, MF, CC, or ALL; `ID`: GO or KEGG IDs; `Description`: description of functional category; `pvalue`: raw p-value of enrichment test; `p.adjust`: p-value adjusted for multiple hypothesis testing based on method specified under `pAdjustMethod`; `qvalue`: q value calculated with R’s `qvalue` function to control FDR; `itemID`: IDs of items (genes for TSEA, drugs for DSEA) overlapping among test and annotation sets; `setSize`: size of the functional category; `GeneRatio`: ratio of genes in the test set that are annotated at a specific GO node or KEGG pathway; `BgRatio`: ratio of background genes that are annotated at a specific GO node or KEGG pathway.

#### 4.1.1.2 With KEGG

The same enrichment test can be performed for KEGG pathways as follows.

    dup_hyperG_k_res <- tsea_dup_hyperG(drugs = drugs, universe = "Default",
                                        type = "KEGG", pvalueCutoff=0.5,
                                        pAdjustMethod="BH", qvalueCutoff = 0.5,
                                        minGSSize = 10, maxGSSize = 500)
    result(dup_hyperG_k_res)

    ## # A tibble: 28 x 9
    ##    ID     Description GeneRatio BgRatio   pvalue p.adjust   qvalue itemID  Count
    ##    <chr>  <chr>       <chr>     <chr>      <dbl>    <dbl>    <dbl> <chr>   <int>
    ##  1 hsa05~ Alcoholism  27/41     187/80~ 3.18e-35 1.18e-33 5.68e-34 55869/~    27
    ##  2 hsa05~ Viral carc~ 26/41     204/80~ 3.04e-32 5.62e-31 2.72e-31 55869/~    26
    ##  3 hsa05~ Huntington~ 17/41     271/80~ 4.25e-15 5.24e-14 2.53e-14 7846/1~    17
    ##  4 hsa04~ Gap juncti~ 12/41     88/8031 8.43e-15 7.80e-14 3.77e-14 7846/1~    12
    ##  5 hsa04~ Phagosome   11/41     152/80~ 1.51e-10 1.11e- 9 5.39e-10 7846/1~    11
    ##  6 hsa05~ Pathogenic~ 11/41     202/80~ 3.17e- 9 1.96e- 8 9.46e- 9 7846/1~    11
    ##  7 hsa04~ Thyroid ho~ 9/41      119/80~ 5.97e- 9 3.16e- 8 1.53e- 8 3065/3~     9
    ##  8 hsa05~ Amphetamin~ 7/41      69/8031 4.52e- 8 2.09e- 7 1.01e- 7 3065/3~     7
    ##  9 hsa04~ Notch sign~ 6/41      53/8031 2.33e- 7 9.57e- 7 4.63e- 7 3065/3~     6
    ## 10 hsa05~ Alzheimer ~ 11/41     334/80~ 5.63e- 7 2.03e- 6 9.80e- 7 7846/1~    11
    ## # ... with 18 more rows

The content of the columns extracted from the `feaResult` object is described under section 4.1.1.1.

### 4.1.2 mGSEA Method

The original GSEA method proposed by Subramanian et al. (2005) uses predefined gene sets \(S\) defined by functional annotation systems such as GO and KEGG. The goal is to determine whether the genes in \(S\) are randomly distributed throughout a ranked test gene list \(L\) (e.g. all genes ranked by LFC) or enriched at the top or bottom of the test list. This is expressed by an Enrichment Score (\(ES\)) reflecting the degree to which a set \(S\) is overrepresented at the extremes of \(L\).

For TSEA, the query is a target set where duplicated entries need to be maintained. To perform GSEA with duplication support, here referred to as mGSEA, the target set is transformed to a score ranked target list \(L_{tar}\) of all targets provided by the corresponding annotation system. For each target in the query target set, its frequency is divided by the number of targets in the target set, which is the weight of that target. For targets present in the annotation system but absent in the target test set, their scores are set to 0. Thus, every target in the annotation system will be assigned a score and then sorted decreasingly to obtain \(L_{tar}\).

In case of TSEA, the original GSEA method cannot be used directly since a large portion of zeros exists in \(L_{tar}\). If the scores of the genes in set \(S\) are all zeros, \(N_R\) (sum of scores of genes in set \(S\)) will be zero, which cannot be used as the denominator. In this case, \(ES\) is set to -1. If only some genes in set \(S\) have scores of zeros then \(N_R\) is set to a larger number to decrease the weight of the genes in \(S\) that have non-zero scores.

\[\begin{equation} N_R=\sum_{g_j\in S}|r_j|^p+min(r_j | r_j > 0)*\sum_{g_j\in S}I_{r_j=0} \end{equation}\]

\(r_j\): score of gene \(j\) in \(L_{tar}\); \(p=1\)

The reason for this modification is that if only one gene in gene set \(S\) has a non-zero score and this gene ranks high in \(L_{tar}\), the weight of this gene will be 1 resulting in an \(ES(S)\) close to 1. Thus, the original GSEA method will score the gene set \(S\) as significantly enriched. However, this is undesirable because in this example only one gene is shared among the target set and the gene set \(S\). Therefore, giving small weights to genes in \(S\) that have scores of zero would decrease the weight of the genes in \(S\) that have scores other than zero, thereby decreasing the false positive rate. To favor truly enriched GO terms and KEGG pathways (gene set \(S\)) at the top of \(L_{tar}\), only gene sets with positive \(ES\) are selected.

#### 4.1.2.1 With GO

The following performs TSEA with the mGSEA method using the same drug test set as in the above `tsea_dup_hyperG` function call. The arguments of the `tsea_mGSEA` function are explained in its help file that can be opened from R with `?tsea_mGSEA`.

    mgsea_res <- tsea_mGSEA(drugs=drugs, type="GO", ont="MF", exponent=1,
                            nPerm=1000, pvalueCutoff=1, minGSSize=5)
    result(mgsea_res)

    ## # A tibble: 93 x 11
    ##    ont   ID    Description setSize enrichmentScore   NES pvalue p.adjust qvalues
    ##    <chr> <chr> <chr>         <int>           <dbl> <dbl>  <dbl>    <dbl>   <dbl>
    ##  1 MF    GO:0~ histone de~      11           1.00  3.81   0.002    0.742   0.742
    ##  2 MF    GO:0~ NAD-depend~      11           1.00  3.81   0.002    0.742   0.742
    ##  3 MF    GO:0~ NAD-depend~      16           0.839 3.72   0.025    0.742   0.742
    ##  4 MF    GO:0~ NAD-depend~      17           0.812 3.69   0.025    0.742   0.742
    ##  5 MF    GO:0~ histone de~      44           0.441 2.52   0.055    0.742   0.742
    ##  6 MF    GO:0~ protein de~      45           0.433 2.50   0.055    0.742   0.742
    ##  7 MF    GO:0~ deacetylas~      58           0.356 2.06   0.07     0.742   0.742
    ##  8 MF    GO:0~ NF-kappaB ~      28           0.265 1.41   0.053    0.742   0.742
    ##  9 MF    GO:0~ hydrolase ~      89           0.250 1.31   0.11     0.742   0.742
    ## 10 MF    GO:0~ potassium ~      10           0.182 0.664  0.445    0.742   0.742
    ## # ... with 83 more rows, and 2 more variables: leadingEdge <chr>,
    ## #   ledge_rank <chr>

The content of the columns extracted from the `feaResult` object is described under section 4.1.1.1. The additional columns specific to the GSEA algorithm are described here.

`enrichmentScore`: \(ES\) from the GSEA algorithm (Subramanian et al. 2005). The score is calculated by walking down the gene list \(L\), increasing a running-sum statistic when we encounter a gene in \(S\) and decreasing when it is not. The magnitude of the increment depends on the gene scores. The \(ES\) is the maximum deviation from zero encountered in the random walk. It corresponds to a weighted Kolmogorov-Smirnov-like statistic.

`NES`: Normalized enrichment score. The positive and negative enrichment scores are normalized separately by permutating the composition of the gene list \(L\) `nPerm` times, and dividing the enrichment score by the mean of the permutation \(ES\) with the same sign.

`pvalue`: The nominal p-value of the \(ES\) is calculated using a permutation test. Specifically, the composition of the gene list \(L\) is permuted and the \(ES\) of the gene set is recomputed for the permutated data generating a null distribution for the ES. The p-value of the observed \(ES\) is then calculated relative to this null distribution.

`leadingEdge`: Genes in the gene set S (functional category) that appear in the ranked list \(L\) at, or before, the point where the running sum reaches its maximum deviation from zero. It can be interpreted as the core of a gene set that accounts for the enrichment signal.

`ledge_rank`: Ranks of genes in ‘leadingEdge’ in gene list \(L\).

#### 4.1.2.2 With KEGG

The same enrichment test can be performed for KEGG pathways as follows.

    mgsea_k_res <- tsea_mGSEA(drugs=drugs, type="KEGG", exponent=1,
                              nPerm=1000, pvalueCutoff=1, minGSSize=2)
    result(mgsea_k_res)

    ## # A tibble: 37 x 10
    ##    ID    Description setSize enrichmentScore   NES pvalue p.adjust qvalues
    ##    <chr> <chr>         <int>           <dbl> <dbl>  <dbl>    <dbl>   <dbl>
    ##  1 hsa0~ Gap juncti~      88          0.135  0.435  0.268    0.799   0.799
    ##  2 hsa0~ Alcoholism      187          0.132  0.287  0.425    0.799   0.799
    ##  3 hsa0~ Viral carc~     204          0.117  0.242  0.45     0.799   0.799
    ##  4 hsa0~ Notch sign~      53          0.105  0.418  0.399    0.799   0.799
    ##  5 hsa0~ Amphetamin~      69          0.0941 0.330  0.422    0.799   0.799
    ##  6 hsa0~ Longevity ~      62          0.0904 0.327  0.434    0.799   0.799
    ##  7 hsa0~ Chronic my~      76          0.0745 0.253  0.508    0.799   0.799
    ##  8 hsa0~ Thyroid ho~     119          0.0716 0.202  0.487    0.799   0.799
    ##  9 hsa0~ Phagosome       152          0.0705 0.174  0.488    0.799   0.799
    ## 10 hsa0~ Huntington~     271          0.0601 0.103  0.61     0.799   0.799
    ## # ... with 27 more rows, and 2 more variables: leadingEdge <chr>,
    ## #   ledge_rank <chr>

The content of the columns extracted from the `feaResult` object is described under sections 4.1.1.1 and 4.1.2.1.

### 4.1.3 MeanAbs Method

The input for the MeanAbs method is \(L_{tar}\), the same as for mGSEA. In this enrichment statistic, \(mabs(S)\), of a gene set \(S\) is calculated as mean absolute scores of the genes in \(S\) (Fang, Tian, and Ji 2012). In order to adjust for size variations in gene set \(S\), 1000 random permutations of \(L_{tar}\) are performed to determine \(mabs(S,\pi)\). Subsequently, \(mabs(S)\) is normalized by subtracting the median of the \(mabs(S,\pi)\) and then dividing by the standard deviation of \(mabs(S,\pi)\) yielding the normalized scores \(Nmabs(S)\). Finally, the portion of \(mabs(S,\pi)\) that is greater than \(mabs(S)\) is used as nominal p-value. The resulting nominal p-values are adjusted for multiple hypothesis testing using the Benjamini-Hochberg method (Benjamini and Hochberg 1995).

#### 4.1.3.1 With GO

The following performs TSEA with the mabs method using the same drug test set as the examples given under the dup_hyperG and tsea_mGSEA sections. The arguments of the `tsea_mabs` function are explained in its help file that can be opened from R with `?tsea_mabs`.

    mabs_res <- tsea_mabs(drugs=drugs, type="GO", ont="MF", nPerm=1000,
                          pvalueCutoff=0.05, minGSSize=5)
    result(mabs_res)

    ## # A tibble: 84 x 10
    ##    ont   ID     Description setSize    mabs Nmabs pvalue p.adjust qvalues itemID
    ##    <chr> <chr>  <chr>         <int>   <dbl> <dbl>  <dbl>    <dbl>   <dbl> <chr>
    ##  1 MF    GO:00~ histone de~      11 0.0576  31.0   0      0       0       HDAC1~
    ##  2 MF    GO:00~ NAD-depend~      11 0.0576  31.0   0      0       0       HDAC1~
    ##  3 MF    GO:00~ NAD-depend~      16 0.0396  30.7   0      0       0       HDAC1~
    ##  4 MF    GO:00~ NAD-depend~      17 0.0373  30.7   0      0       0       HDAC1~
    ##  5 MF    GO:00~ histone de~      44 0.0144  28.4   0      0       0       HDAC1~
    ##  6 MF    GO:00~ protein de~      45 0.0141  28.4   0      0       0       HDAC1~
    ##  7 MF    GO:00~ deacetylas~      58 0.0109  27.3   0      0       0       HDAC1~
    ##  8 MF    GO:00~ NF-kappaB ~      28 0.00784 23.4   0      0       0       HDAC1~
    ##  9 MF    GO:00~ hydrolase ~      89 0.00713 24.8   0      0       0       HDAC1~
    ## 10 MF    GO:00~ potassium ~      10 0.00488  8.01  0.003  0.00951 0.00191 HDAC4
    ## # ... with 74 more rows

The content of the columns extracted from the `feaResult` object is explained under section 4.1.1.1. The columns specific to the mabs algorithm are described below.

`mabs`: Given a scored ranked gene list \(L\), \(mabs(S)\) represents the mean absolute score of the genes in set \(S\).

`Nmabs`: \(mabs(S)\) normalized.

#### 4.1.3.2 With KEGG

The same enrichment test can be performed for KEGG pathways as follows.

    mabs_k_res <- tsea_mabs(drugs=drugs, type="KEGG", nPerm=1000,
                            pvalueCutoff=0.2, minGSSize=5)
    result(mabs_k_res)

    ## # A tibble: 31 x 9
    ##    ID     Description   setSize    mabs Nmabs pvalue p.adjust qvalues itemID
    ##    <chr>  <chr>           <int>   <dbl> <dbl>  <dbl>    <dbl>   <dbl> <chr>
    ##  1 hsa05~ Alcoholism        187 0.00352 16.6       0        0       0 10013/100~
    ##  2 hsa04~ Gap junction       88 0.00333 12.1       0        0       0 10376/103~
    ##  3 hsa05~ Viral carcin~     204 0.00311 16.0       0        0       0 10013/100~
    ##  4 hsa04~ Notch signal~      53 0.00276  8.61      0        0       0 3065/3066
    ##  5 hsa05~ Amphetamine ~      69 0.00247  8.01      0        0       0 1812/3065~
    ##  6 hsa04~ Longevity re~      62 0.00236  7.62      0        0       0 3065/3066
    ##  7 hsa05~ Chronic myel~      76 0.00193  7.20      0        0       0 3065/3066
    ##  8 hsa04~ Thyroid horm~     119 0.00184  8.66      0        0       0 3065/3066~
    ##  9 hsa04~ Phagosome         152 0.00177  8.74      0        0       0 10376/103~
    ## 10 hsa05~ Huntington d~     271 0.00153  8.80      0        0       0 10376/103~
    ## # ... with 21 more rows

## 4.2 DSEA

Instead of translating ranked lists of drugs into target sets, as for TSEA, the functional annotation categories of the targets can be assigned to the drugs directly to perform Drug Set Enrichment Analysis (DSEA) instead. Since the
drug lists from GESS results are usually unique, this strategy overcomes the duplication problem of the TSEA approach. This way classical enrichment methods, such as GSEA or tests based on the hypergeometric distribution, can be readily applied without major modifications to the underlying statistical methods. As explained above, TSEA and DSEA performed with the same enrichment statistics are not expected to generate identical results. Rather they often complement each other’s strengths and weaknesses.

### 4.2.1 Hypergeometric Test

The following performs DSEA with `signatureSearch's` hypergeometric test function called `dsea_hyperG` using the same drug test set as the examples given under in the TSEA section. The arguments are explained in its help file that can be opened from R with `?dsea_hyperG`.

#### 4.2.1.1 With GO

As functional annotation system the following DSEA example uses GO.

    drugs <- unique(result(lincs)$pert[1:10])
    hyperG_res <- dsea_hyperG(drugs=drugs, type="GO", ont="MF")
    result(hyperG_res)

    ## # A tibble: 51 x 10
    ##    ont   ID    Description GeneRatio BgRatio   pvalue p.adjust  qvalue itemID
    ##    <chr> <chr> <chr>       <chr>     <chr>      <dbl>    <dbl>   <dbl> <chr>
    ##  1 MF    GO:0~ Krueppel-a~ 4/8       27/225~ 1.13e-10  5.89e-8 1.73e-8 vorin~
    ##  2 MF    GO:0~ histone de~ 5/8       199/22~ 2.78e- 9  3.04e-7 8.93e-8 vorin~
    ##  3 MF    GO:0~ NAD-depend~ 5/8       199/22~ 2.78e- 9  3.04e-7 8.93e-8 vorin~
    ##  4 MF    GO:0~ NAD-depend~ 5/8       261/22~ 1.08e- 8  3.38e-7 9.93e-8 vorin~
    ##  5 MF    GO:0~ NAD-depend~ 5/8       266/22~ 1.19e- 8  3.38e-7 9.93e-8 vorin~
    ##  6 MF    GO:0~ nucleosoma~ 4/8       88/225~ 1.49e- 8  3.88e-7 1.14e-7 vorin~
    ##  7 MF    GO:0~ core promo~ 4/8       94/225~ 1.95e- 8  4.61e-7 1.35e-7 vorin~
    ##  8 MF    GO:0~ histone de~ 5/8       318/22~ 2.91e- 8  6.11e-7 1.79e-7 vorin~
    ##  9 MF    GO:0~ protein de~ 5/8       323/22~ 3.15e- 8  6.11e-7 1.79e-7 vorin~
    ## 10 MF    GO:0~ nucleosome~ 4/8       106/22~ 3.17e- 8  6.11e-7 1.79e-7 vorin~
    ## # ... with 41 more rows, and 1 more variable: Count <int>

#### 4.2.1.2 With KEGG

The same DSEA test can be performed for KEGG pathways as follows.

    hyperG_k_res <- dsea_hyperG(drugs = drugs, type = "KEGG",
                                pvalueCutoff = 1, qvalueCutoff = 1,
                                minGSSize = 10, maxGSSize = 2000)
    result(hyperG_k_res)

    ## # A tibble: 34 x 9
    ##    ID     Description  GeneRatio BgRatio  pvalue p.adjust  qvalue itemID   Count
    ##    <chr>  <chr>        <chr>     <chr>     <dbl>    <dbl>   <dbl> <chr>    <int>
    ##  1 hsa04~ Notch signa~ 4/8       146/20~ 1.70e-7  5.79e-6 3.76e-6 vorinos~     4
    ##  2 hsa05~ Alcoholism   6/8       1131/2~ 7.15e-7  1.22e-5 7.91e-6 vorinos~     6
    ##  3 hsa05~ Amphetamine~ 5/8       916/20~ 8.89e-6  1.01e-4 6.55e-5 vorinos~     5
    ##  4 hsa05~ Huntington ~ 5/8       1085/2~ 2.03e-5  1.73e-4 1.12e-4 vorinos~     5
    ##  5 hsa05~ Viral carci~ 5/8       1175/2~ 2.99e-5  2.04e-4 1.32e-4 vorinos~     5
    ##  6 hsa04~ Longevity r~ 4/8       578/20~ 4.03e-5  2.28e-4 1.48e-4 vorinos~     4
    ##  7 hsa04~ Cell cycle   4/8       730/20~ 1.00e-4  4.87e-4 3.16e-4 vorinos~     4
    ##  8 hsa05~ Chronic mye~ 4/8       822/20~ 1.59e-4  6.75e-4 4.39e-4 vorinos~     4
    ##  9 hsa05~ Transcripti~ 4/8       996/20~ 3.33e-4  1.26e-3 8.19e-4 vorinos~     4
    ## 10 hsa04~ Thyroid hor~ 4/8       1175/2~ 6.28e-4  2.14e-3 1.39e-3 vorinos~     4
    ## # ... with 24 more rows

The content of the columns extracted from the feaResult object is described under section 4.1.1.1.

### 4.2.2 GSEA Method

The following performs DSEA with the GSEA method using as test set drug labels combined with scores. Instead of using only the drug labels in the test set, the GSEA method requires the labels as well as the scores used for ranking the drug list in the GESS result. The scores are usually the similarity metric used to rank the results of the corresponding GESS method, here the NCS values from the LINCS method. The arguments of the `dsea_GSEA` function are explained in its help file that can be opened from R with `?dsea_GSEA`.

#### 4.2.2.1 With GO

As functional annotation system the following DSEA example uses GO.

    dl <- abs(result(lincs)$NCS); names(dl) <- result(lincs)$pert
    dl <- dl[dl>0]
    dl <- dl[!duplicated(names(dl))]
    gsea_res <- dsea_GSEA(drugList=dl, type="GO", ont="MF", exponent=1, nPerm=1000,
                          pvalueCutoff=0.2, minGSSize=5)

    ## Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize
    ## = minGSSize, : You are trying to run fgseaSimple. It is recommended to use
    ## fgseaMultilevel. To run fgseaMultilevel, you need to remove the nperm argument
    ## in the fgsea function call.

    ## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
    ## gseaParam, : All values in the stats vector are greater than zero and scoreType
    ## is "std", maybe you should switch to scoreType = "pos".

    result(gsea_res)

    ## # A tibble: 55 x 11
    ##    ont   ID    Description setSize enrichmentScore   NES  pvalue p.adjust
    ##    <chr> <chr> <chr>         <int>           <dbl> <dbl>   <dbl>    <dbl>
    ##  1 MF    GO:0~ histone de~       5           0.951  2.24 0.00149   0.0254
    ##  2 MF    GO:0~ hydrolase ~       5           0.951  2.24 0.00149   0.0254
    ##  3 MF    GO:0~ NAD-depend~       5           0.951  2.24 0.00149   0.0254
    ##  4 MF    GO:0~ deacetylas~       5           0.951  2.24 0.00149   0.0254
    ##  5 MF    GO:0~ histone de~       5           0.951  2.24 0.00149   0.0254
    ##  6 MF    GO:0~ NAD-depend~       5           0.951  2.24 0.00149   0.0254
    ##  7 MF    GO:0~ protein de~       5           0.951  2.24 0.00149   0.0254
    ##  8 MF    GO:0~ NAD-depend~       5           0.951  2.24 0.00149   0.0254
    ##  9 MF    GO:0~ transcript~       5           0.931  2.20 0.00149   0.0254
    ## 10 MF    GO:0~ hydrolase ~       6           0.863  2.17 0.00143   0.0254
    ## # ... with 45 more rows, and 3 more variables: qvalues <dbl>,
    ## #   leadingEdge <chr>, ledge_rank <chr>

#### 4.2.2.2 With KEGG

The same DSEA test can be performed for KEGG pathways as follows.

    gsea_k_res <- dsea_GSEA(drugList=dl, type="KEGG", exponent=1, nPerm=1000,
                            pvalueCutoff=1, minGSSize=5)

    ## Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize
    ## = minGSSize, : You are trying to run fgseaSimple. It is recommended to use
    ## fgseaMultilevel. To run fgseaMultilevel, you need to remove the nperm argument
    ## in the fgsea function call.

    ## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
    ## gseaParam, : All values in the stats vector are greater than zero and scoreType
    ## is "std", maybe you should switch to scoreType = "pos".

    result(gsea_k_res)

    ## # A tibble: 62 x 10
    ##    ID    Description setSize enrichmentScore   NES  pvalue p.adjust qvalues
    ##    <chr> <chr>         <int>           <dbl> <dbl>   <dbl>    <dbl>   <dbl>
    ##  1 hsa0~ Amphetamin~       8           0.785  2.18 0.00134   0.0363  0.0302
    ##  2 hsa0~ Alcoholism       10           0.731  2.16 0.00132   0.0363  0.0302
    ##  3 hsa0~ Cell cycle        6           0.838  2.09 0.00288   0.0363  0.0302
    ##  4 hsa0~ Epstein-Ba~       6           0.838  2.09 0.00288   0.0363  0.0302
    ##  5 hsa0~ Transcript~       7           0.772  2.04 0.00410   0.0363  0.0302
    ##  6 hsa0~ Longevity ~       7           0.762  2.02 0.00410   0.0363  0.0302
    ##  7 hsa0~ Huntington~       8           0.718  1.99 0.00401   0.0363  0.0302
    ##  8 hsa0~ Chronic my~       8           0.717  1.99 0.00535   0.0368  0.0306
    ##  9 hsa0~ Thyroid ho~       9           0.673  1.93 0.00528   0.0368  0.0306
    ## 10 hsa0~ Viral carc~      10           0.639  1.89 0.00789   0.0489  0.0407
    ## # ... with 52 more rows, and 2 more variables: leadingEdge <chr>,
    ## #   ledge_rank <chr>

Since the annotation system are drug-to-functional category mappings, the “leadingEdge” column contains identifiers of drugs instead of targets.

## 4.3 Comparing FEA Results

The `comp_fea_res` function re-ranks the functional categories across the different FEA methods by using the mean rank of each functional category across the 5 FEA methods. Here the functional categories are re-ranked by their mean rank values in increasing order. Since the functional categories are not always present in all enrichment results, the mean rank of a functional category is corrected by an adjustment factor that is the number of enrichment result methods used divided by the number of occurrences of a functional category. For instance, if a functional category is only present in the result of one method, its mean rank will be increased accordingly. The following plots use the `pvalue` column in the result tables for this ranking approach. Alternative columns can be chosen under the `rank_stat` argument. After re-ranking only the top ranking functional categories are shown, here 20. Their number can be changed under the `Nshow` argument.

    table_list = list("dup_hyperG" = result(dup_hyperG_res),
                      "mGSEA" = result(mgsea_res),
                      "mabs" = result(mabs_res),
                      "hyperG" = result(hyperG_res),
                      "GSEA" = result(gsea_res))
    comp_fea_res(table_list, rank_stat="pvalue", Nshow=20)

Figure 4: Ranking comparison of top GO terms across FEA methods. The dots in the plot represent the p-values of the 20 top ranking MF GO terms (y-axis) with heat color coding according to the color gradient in the legend on the right. The GO terms have been ordered from the top to the bottom of the plot by increasing mean rank values calculated for each GO term across the 5 FEA methods (x-axis).

    table_list = list("dup_hyperG" = result(dup_hyperG_k_res),
                      "mGSEA" = result(mgsea_k_res),
                      "mabs" = result(mabs_k_res),
                      "hyperG" = result(hyperG_k_res),
                      "GSEA" = result(gsea_k_res))
    comp_fea_res(table_list, rank_stat="pvalue", Nshow=20)

Figure 5: Ranking comparison among top KEGG pathways across FEA methods. The plot depicts the same information as the previous plot but using KEGG pathways instead of GO terms as functional annotation system.

## 4.4 Summary of FEA Results

The enrichment rankings of the functional categories (GO and KEGG) show a reasonable degree of agreement among the five FEA methods. For instance, all five showed a high degree of enrichment for histone deacetylase pathways that are indeed targeted by one of the query drugs, here vorinostat. Since each method has its strengths and weaknesses, the usage of a consensus approach could be considered by combining the rankings of functional categories from all or several FEA methods.

# 5 Drug-Target Network Visualization

Functional modules of the above GESS and FEA results can be rendered as interactive drug-target networks using the `dtnetplot` function from `signatureSearch`. For this, a character vector of drug names along with an identifier of a chosen functional category are passed on to the `drugs` and `set` arguments, respectively. The resulting plot depicts the corresponding drug-target interaction network. Its interactive features allow the user to zoom in and out of the network, and to select network components in the drop-down menu located in the upper left corner of the plot.

## 5.1 With GO

The following example demonstrates how to construct a drug-target network for two GO categories that are enriched in the results obtained in the FEA section of this vignette. This way one can visualize drug sets in the context of the cellular networks and pathways they affect.

Network for NAD-dependent histone deacetylase activity (GO:0032041)

    dtnetplot(drugs = drugs(dup_hyperG_res), set = "GO:0032041", ont = "MF",
              desc="NAD-dependent histone deacetylase activity (H3-K14 specific)")

Figure 6: Drug-target network for NAD-dependent histone deacetylase activity. The given network graph illustrates the interactions among drugs and their target proteins in the chosen pathway. Only drugs with at least one target protein are included. The nodes in the network represent drugs and targets depicted as squares and circles, respectively. The interactions among drugs and targets are indicated by non-weighted lines (yes/no relationship). The color of the nodes is defined by the number of their connections to other nodes.

Network for NF-kappaB binding (GO:0051059)

    dtnetplot(drugs = drugs(dup_hyperG_res), set = "GO:0051059", ont = "MF",
              desc="NF-kappaB binding")

Figure 7: Drug-target network for NF-kappaB binding. The details of the plot are given in the legend of Figure 6.

## 5.2 With KEGG

The same drug-target network plots can be rendered for KEGG pathways as follows.

Network for KEGG pathway: Alcoholism (hsa05034)

    dtnetplot(drugs = drugs(dup_hyperG_k_res), set = "hsa05034",
              desc="Alcoholism")

Figure 8: Drug-target network for Alcoholism. The details of the plot are given in the legend of Figure 6.

Network for KEGG pathway: Longevity regulating pathway (hsa04213)

    dtnetplot(drugs = drugs, set = "hsa04213",
              desc="Longevity regulating pathway - multiple species")

Figure 9: Drug-target network for Longevity regulating pathway. The details of the plot are given in the legend of Figure 6.

# 6 Run Workflow

The `runWF` function could run the entire GESS/FEA workflow automatically when providing the query drug and cell type, as well as selecting the reference database (e.g. `lincs` or path to the custom reference database), defining the specific GESS and FEA methods. When the query drug and cell type were provided, the query GES was internally drawn from the reference database. The N (defined by the `N_gess_drugs` argument) top ranking hits in the GESS tables were then used for FEA where three different annotation systems were used: GO Molecular Function (GO MF), GO Biological Process (GO BP) and KEGG pathways.

The GESS/FEA results will be stored in a list object in R session. A working environment named by the use case will be created under users current working directory or under other directory defined by users. This environment contains a folder where the GESS/FEA result tables were written to. The working environment also contains a template Rmd vignette as well as a rended HTML report, users could make modifications on the Rmd vignette as they need and re-render it to generate their HTML report by running `rmarkdown::render("GESS_FEA_report.Rmd")` in R session or `Rscript -e "rmarkdown::render('GESS_FEA_report.Rmd')"` from bash commandline.

    drug <- "vorinostat"; cell <- "SKB"
    refdb <- system.file("extdata", "sample_db.h5", package="signatureSearch")
    env_dir <- tempdir()
    wf_list <- runWF(drug, cell, refdb, gess_method="LINCS",
        fea_method="dup_hyperG", N_gess_drugs=10, env_dir=env_dir, tau=FALSE)

# 7 Supplemental Material

## 7.1 Construction of Toy Database

Access to the LINCS database is provided via the associated signatureSearchData package hosted on Bioconductor’s `ExperimentHub`. The following provides the code for constructing the toy database used by the sample code of this vignette. To save time building this vignette, the evaluated components of its sample code use a pre-generated instance of the toy database that is stored in the `extdata` directory of the `signatureSearch` package. Thus, the following code section is not evaluated. It also serves as an example how to construct other custom instances of the LINCS database. Additional details on this topic are provided in the vignette of the `signatureSearchData` package.

    library(rhdf5)
    eh <- ExperimentHub::ExperimentHub()
    lincs <- eh[["EH3226"]]
    hdacs <- c("vorinostat","trichostatin-a","pyroxamide","HC-toxin","rhamnetin")
    hdacs_trts <- paste(hdacs, "SKB", "trt_cp", sep="__")
    all_trts <- drop(h5read(lincs, "colnames"))
    # Select treatments in SKB cell and not BRD compounds
    all_trts2 <- all_trts[!grepl("BRD-", all_trts) & grepl("__SKB__", all_trts)]
    set.seed(11)
    rand_trts <- sample(setdiff(all_trts2, hdacs_trts), 95)
    toy_trts <- c(hdacs_trts, rand_trts)
    library(SummarizedExperiment); library(HDF5Array)
    toy_db <- SummarizedExperiment(HDF5Array(lincs, name="assay"))
    rownames(toy_db) <- HDF5Array(db_path, name="rownames")
    colnames(toy_db) <- HDF5Array(db_path, name="colnames")
    toy_db <- round(as.matrix(assay(toy_db)[,toy_trts]),2)
    set.seed(11)
    gene_idx <- sample.int(nrow(toy_db),5000)
    toy_db2 <- toy_db[gene_idx,]
    # The sample_db is stored in the current directory of user's R session
    getwd()
    createEmptyH5("sample_db.h5", level=9, delete_existing=TRUE)
    append2H5(toy_db2, "sample_db.h5")
    h5ls("sample_db.h5")

# 8 Session Info

    sessionInfo()

    ## R version 4.0.0 beta (2020-04-14 r78227)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] org.Hs.eg.db_3.10.0         AnnotationDbi_1.49.1
    ##  [3] signatureSearchData_1.1.0   HDF5Array_1.15.18
    ##  [5] rhdf5_2.31.10               ggplot2_3.3.0
    ##  [7] signatureSearch_1.1.5       SummarizedExperiment_1.17.5
    ##  [9] DelayedArray_0.13.12        matrixStats_0.56.0
    ## [11] Biobase_2.47.3              GenomicRanges_1.39.3
    ## [13] GenomeInfoDb_1.23.17        IRanges_2.21.8
    ## [15] S4Vectors_0.25.15           BiocGenerics_0.33.3
    ## [17] Rcpp_1.0.4.6
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] AnnotationHub_2.19.11         fastmatch_1.1-0
    ##   [3] BiocFileCache_1.11.6          plyr_1.8.6
    ##   [5] igraph_1.2.5                  GSEABase_1.49.1
    ##   [7] splines_4.0.0                 BiocParallel_1.21.2
    ##   [9] urltools_1.7.3                digest_0.6.25
    ##  [11] htmltools_0.4.0               GOSemSim_2.13.1
    ##  [13] viridis_0.5.1                 GO.db_3.10.0
    ##  [15] fansi_0.4.1                   magrittr_1.5
    ##  [17] memoise_1.1.0                 limma_3.43.8
    ##  [19] readr_1.3.1                   annotate_1.65.1
    ##  [21] graphlayouts_0.6.0            R.utils_2.9.2
    ##  [23] enrichplot_1.7.4              prettyunits_1.1.1
    ##  [25] colorspace_1.4-1              blob_1.2.1
    ##  [27] rappdirs_0.3.1                ggrepel_0.8.2
    ##  [29] xfun_0.13                     dplyr_0.8.5
    ##  [31] crayon_1.3.4                  RCurl_1.98-1.1
    ##  [33] jsonlite_1.6.1                graph_1.65.3
    ##  [35] scatterpie_0.1.4              glue_1.4.0
    ##  [37] polyclip_1.10-0               gtable_0.3.0
    ##  [39] zlibbioc_1.33.1               XVector_0.27.2
    ##  [41] Rhdf5lib_1.9.3                scales_1.1.0
    ##  [43] DOSE_3.13.2                   DBI_1.1.0
    ##  [45] viridisLite_0.3.0             xtable_1.8-4
    ##  [47] progress_1.2.2                gridGraphics_0.5-0
    ##  [49] bit_1.1-15.2                  europepmc_0.3
    ##  [51] preprocessCore_1.49.2         htmlwidgets_1.5.1
    ##  [53] httr_1.4.1                    fgsea_1.13.5
    ##  [55] RColorBrewer_1.1-2            ellipsis_0.3.0
    ##  [57] pkgconfig_2.0.3               XML_3.99-0.3
    ##  [59] R.methodsS3_1.8.0             farver_2.0.3
    ##  [61] dbplyr_1.4.2                  utf8_1.1.4
    ##  [63] labeling_0.3                  ggplotify_0.0.5
    ##  [65] tidyselect_1.0.0              rlang_0.4.5
    ##  [67] reshape2_1.4.4                later_1.0.0
    ##  [69] munsell_0.5.0                 BiocVersion_3.11.1
    ##  [71] tools_4.0.0                   visNetwork_2.0.9
    ##  [73] downloader_0.4                cli_2.0.2
    ##  [75] RSQLite_2.2.0                 ExperimentHub_1.13.7
    ##  [77] ggridges_0.5.2                evaluate_0.14
    ##  [79] stringr_1.4.0                 fastmap_1.0.1
    ##  [81] yaml_2.2.1                    knitr_1.28
    ##  [83] bit64_0.9-7                   tidygraph_1.1.2
    ##  [85] purrr_0.3.3                   ggraph_2.0.2
    ##  [87] mime_0.9                      R.oo_1.23.0
    ##  [89] DO.db_2.9                     xml2_1.3.1
    ##  [91] compiler_4.0.0                curl_4.3
    ##  [93] interactiveDisplayBase_1.25.0 affyio_1.57.0
    ##  [95] tibble_3.0.0                  tweenr_1.0.1
    ##  [97] stringi_1.4.6                 lattice_0.20-41
    ##  [99] Matrix_1.2-18                 vctrs_0.2.4
    ## [101] pillar_1.4.3                  lifecycle_0.2.0
    ## [103] BiocManager_1.30.10           triebeard_0.3.0
    ## [105] data.table_1.12.8             cowplot_1.0.0
    ## [107] bitops_1.0-6                  httpuv_1.5.2
    ## [109] qvalue_2.19.0                 R6_2.4.1
    ## [111] affy_1.65.1                   promises_1.1.0
    ## [113] gridExtra_2.3                 MASS_7.3-51.5
    ## [115] assertthat_0.2.1              withr_2.1.2
    ## [117] GenomeInfoDbData_1.2.3        hms_0.5.3
    ## [119] clusterProfiler_3.15.4        grid_4.0.0
    ## [121] tidyr_1.0.2                   rmarkdown_2.1
    ## [123] rvcheck_0.1.8                 ggforce_0.3.1
    ## [125] shiny_1.4.0.2

# 9 Funding

This project is funded by NIH grants U19AG02312 and U24AG051129 awarded by the National Institute on Aging (NIA). Subcomponents of the environment are based on methods developed by projects funded by NSF awards ABI-1661152 and PGRP-1810468. The High-Performance Computing (HPC) resources used for testing and optimizing the code of this project were funded by NIH and NSF grants 1S10OD016290-01A1 and MRI-1429826, respectively.

# References

Benjamini, Y, and Y Hochberg. 1995. “Controlling the false discovery rate - a practical and powerful approach to multiple testing.” J. R. Stat. Soc. Series B Stat. Methodol. 57: 289–300.

Edgar, Ron, Michael Domrachev, and Alex E Lash. 2002. “Gene Expression Omnibus: NCBI Gene Expression and Hybridization Array Data Repository.” Nucleic Acids Res. 30 (1): 207–10. https://www.ncbi.nlm.nih.gov/pubmed/11752295.

Fang, Zhaoyuan, Weidong Tian, and Hongbin Ji. 2012. “A Network-Based Gene-Weighting Approach for Pathway Analysis.” Cell Res. 22 (3): 565–80. http://dx.doi.org/10.1038/cr.2011.149.

Gaulton, Anna, Anne Hersey, Michał Nowotka, A Patrı́cia Bento, Jon Chambers, David Mendez, Prudence Mutowo, et al. 2017. “The ChEMBL database in 2017.” Nucleic Acids Res. 45 (D1): D945–D954. https://doi.org/10.1093/nar/gkw1074.

Graham J. G. Upton. 1992. “Fisher’s Exact Test.” J. R. Stat. Soc. Ser. A Stat. Soc. 155 (3): 395–402. http://www.jstor.org/stable/2982890.

Kuhn, Michael, Damian Szklarczyk, Andrea Franceschini, Monica Campillos, Christian von Mering, Lars Juhl Jensen, Andreas Beyer, and Peer Bork. 2010. “STITCH 2: An Interaction Network Database for Small Molecules and Proteins.” Nucleic Acids Res. 38 (Database issue): D552–6. http://dx.doi.org/10.1093/nar/gkp937.

Lamb, Justin, Emily D Crawford, David Peck, Joshua W Modell, Irene C Blat, Matthew J Wrobel, Jim Lerner, et al. 2006. “The Connectivity Map: Using Gene-Expression Signatures to Connect Small Molecules, Genes, and Disease.” Science 313 (5795): 1929–35. http://dx.doi.org/10.1126/science.1132939.

Peck, David, Emily D Crawford, Kenneth N Ross, Kimberly Stegmaier, Todd R Golub, and Justin Lamb. 2006. “A Method for High-Throughput Gene Expression Signature Analysis.” Genome Biol. 7 (7): R61. http://dx.doi.org/10.1186/gb-2006-7-7-r61.

Sandmann, Thomas, Sarah K Kummerfeld, Robert Gentleman, and Richard Bourgon. 2014. “gCMAP: User-Friendly Connectivity Mapping with R.” Bioinformatics 30 (1): 127–28. http://dx.doi.org/10.1093/bioinformatics/btt592.

Subramanian, Aravind, Rajiv Narayan, Steven M Corsello, David D Peck, Ted E Natoli, Xiaodong Lu, Joshua Gould, et al. 2017. “A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles.” Cell 171 (6): 1437–1452.e17. http://dx.doi.org/10.1016/j.cell.2017.10.049.

Subramanian, Aravind, Pablo Tamayo, Vamsi K Mootha, Sayan Mukherjee, Benjamin L Ebert, Michael A Gillette, Amanda Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proc. Natl. Acad. Sci. U. S. A. 102 (43): 15545–50. http://dx.doi.org/10.1073/pnas.0506580102.

Wishart, David S, Yannick D Feunang, An C Guo, Elvis J Lo, Ana Marcu, Jason R Grant, Tanvir Sajed, et al. 2018. “DrugBank 5.0: A Major Update to the DrugBank Database for 2018.” Nucleic Acids Res. 46 (D1): D1074–D1082. http://dx.doi.org/10.1093/nar/gkx1037.