# Gene Set Enrichment Analysis with Networks

#### Dongmin Jung

#### January 24, 2018

# 1 Introduction

Biological molecules in a living organism seldom work individually. They usually interact each other in a cooperative way. Biological process is too complicated to understand without considering such interactions. These associations can be conveniently represented as networks (graphs) with nodes (vertices) for molecules and links (edges) for interactions between them. Thus, network-based procedures can be seen as powerful methods for studying complex process. For example, cancer may be better characterized by frequently mutated or dysregulated pathways than driver mutations. However, many methods are devised for analyzing individual genes. It is said that techniques based on biological networks such as gene co-expression are more precise ways to represent information than those using lists of genes only. This package is aimed to integrate the gene expression and biological network. A biological network is constructed from gene expression data and it is used for Gene Set Enrichment Analysis.

Recently, network module-based methods have been proposed as a powerful tool to analyze and interpret gene expression data. For example, a co-expression gene network exhibit distinctive connection structures, including their topologies. In particular, the highly connected genes, or hub genes located at the functional center of the network, should have a high tendency to be associated with their phenotypic outcome. In case that a pathway is represented as a network with nodes and edges, its topology is essential for evaluating the importance of the pathway. However, the widely used statistical methods were designed for individual genes, which may detect too many irrelevant significantly genes or too few genes to describe the phenotypic changes. One possible reason is that most of the existing functional enrichment analyses ignore topological information. Therfore, they are limited for applicability. To overcome this limitation, the original functional enrichment method is extended by assigning network centralities such as node weights. The rationale behind this is that all nodes are not equally important because of hub genes. A fundamental problem in the analysis of complex networks is the determination of how important a certain node is. Centrality measures such as node strength show how influential that node is in the overall network. It assigns relative scores to all nodes within the graph. For simplicity, degree centrality or node strength is used as a centrality measure by default. For over-representation analysis, the label propagation algorithm is employed to score nodes in this package.

# 2 Example

## 2.1 GSEA with the list of genes, based on the label propagation

Label propagation is a semi-supervised learning that assigns labels for unlabelled nodes in a network. Each node will be labeled, based on the corresponding score. These score can be used as a statistic for GSEA. Thus ORA can be performed by GSEA. For the example of GSEA for the label propagation, consider Lung Adenocarcinoma data from TCGA. We have set our threshold to 0.8, so only those nodes with similarities higher than 0.8 will be considered neighbors. Consider recurrently mutated genes from MutSig and KEGG pathways. In the result, we can see that these mutated genes are enriched at cancer relevant signaling pathways such as ERBB and MAPK pathways.

    library(gsean)
    library(TCGAbiolinks)
    library(DESeq2)
    library(PPInfer)
    # TCGA LUAD
    query <- GDCquery(project = "TCGA-LUAD",
                      data.category = "Gene expression",
                      data.type = "Gene expression quantification",
                      platform = "Illumina HiSeq",
                      file.type  = "normalized_results",
                      experimental.strategy = "RNA-Seq",
                      legacy = TRUE)
    GDCdownload(query, method = "api")
    invisible(capture.output(data <- GDCprepare(query)))
    exprs.LUAD <- assay(data)
    # remove duplicated gene names
    exprs.LUAD <- exprs.LUAD[-which(duplicated(rownames(exprs.LUAD))),]
    # list of genes
    recur.mut.gene <- c("KRAS", "TP53", "STK11", "RBM10", "SPATA31C1", "KRTAP4-11",
                        "DCAF8L2", "AGAP6", "KEAP1", "SETD2", "ZNF679", "FSCB",
                        "BRAF", "ZNF770", "U2AF1", "SMARCA4", "HRNR", "EGFR")

    # KEGG_hsa
    load(system.file("data", "KEGG_hsa.rda", package = "gsean"))

    # GSEA
    set.seed(1)
    result.GSEA <- gsean(KEGG_hsa, recur.mut.gene, exprs.LUAD, threshold = 0.8)
    invisible(capture.output(p <- GSEA.barplot(result.GSEA, category = 'pathway',
                                               score = 'NES', pvalue = 'padj',
                                               sort = 'padj', top = 20)))
    p <- p + scale_fill_gradient(low = "red", high = "blue") +
      theme(plot.margin = margin(10, 10, 10, 75))
    plotly::ggplotly(p + guides(fill = FALSE))

## 2.2 GSEA with statistics, based on the degree centrality

Gene expression data are analyzed to identify differentially expressed genes. Consider pasilla RNA-seq count data. By using the Wald statistic in this example, GSEA can be performed with Gene Ontology terms from http://www.go2msig.org/cgi-bin/prebuilt.cgi?taxid=7227. Thus, it is expected that we may find the biological functions related to change in phenotype from the network, rather than individual genes.

    library(gsean)
    library(pasilla)
    library(DESeq2)
    library(PPInfer)
    # pasilla count data
    pasCts <- system.file("extdata", "pasilla_gene_counts.tsv",
                          package = "pasilla", mustWork = TRUE)
    cts <- as.matrix(read.csv(pasCts, sep="\t", row.names = "gene_id"))
    condition <- factor(c(rep("untreated", 4), rep("treated", 3)))
    dds <- DESeqDataSetFromMatrix(
      countData = cts,
      colData   = data.frame(condition),
      design    = ~ 0 + condition)
    # filtering
    keep <- rowSums(counts(dds)) >= 10
    dds <- dds[keep,]
    # differentially expressed genes
    dds <- DESeq(dds)
    resultsNames(dds)

    ## [1] "conditiontreated"   "conditionuntreated"

    res <- results(dds, contrast = list("conditiontreated", "conditionuntreated"), listValues = c(1, -1))
    statistic <- res$stat
    names(statistic) <- rownames(res)
    exprs.pasilla <- counts(dds, normalized = TRUE)

    # convert gene id
    library(org.Dm.eg.db)
    gene.id <- AnnotationDbi::select(org.Dm.eg.db, names(statistic), "ENTREZID", "FLYBASE")
    names(statistic) <- gene.id[,2]
    rownames(exprs.pasilla) <- gene.id[,2]

    # GO_dme
    load(system.file("data", "GO_dme.rda", package = "gsean"))

    # GSEA
    set.seed(1)
    result.GSEA <- gsean(GO_dme, statistic, exprs.pasilla)

    ## log2 transformed for correlation

    invisible(capture.output(p <- GSEA.barplot(result.GSEA, category = 'pathway',
                                               score = 'NES', top = 50, pvalue = 'padj',
                                               sort = 'padj', decreasing = FALSE,
                                               numChar = 110)))
    p <- p + scale_fill_continuous(low = 'red', high = 'green') +
      theme(plot.margin = margin(10, 10, 10, 50))
    plotly::ggplotly(p)

# 3 Session info

    sessionInfo()

    ## R version 3.5.0 (2018-04-23)
    ## 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] stats4    parallel  stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] org.Dm.eg.db_3.6.0          AnnotationDbi_1.42.0
    ##  [3] DESeq2_1.20.0               SummarizedExperiment_1.10.0
    ##  [5] DelayedArray_0.6.0          BiocParallel_1.14.0
    ##  [7] matrixStats_0.53.1          Biobase_2.40.0
    ##  [9] GenomicRanges_1.32.0        GenomeInfoDb_1.16.0
    ## [11] IRanges_2.14.0              S4Vectors_0.18.0
    ## [13] pasilla_1.7.0               gsean_1.0.0
    ## [15] PPInfer_1.6.0               yeastExpData_0.25.0
    ## [17] graph_1.58.0                BiocGenerics_0.26.0
    ## [19] STRINGdb_1.20.0             igraph_1.2.1
    ## [21] ggplot2_2.2.1               kernlab_0.9-26
    ## [23] biomaRt_2.36.0              fgsea_1.6.0
    ## [25] Rcpp_0.12.16
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] backports_1.1.2        Hmisc_4.1-1            fastmatch_1.1-0
    ##   [4] plyr_1.8.4             lazyeval_0.2.1         splines_3.5.0
    ##   [7] crosstalk_1.0.0        robust_0.4-18          digest_0.6.15
    ##  [10] foreach_1.4.4          htmltools_0.3.6        GO.db_3.6.0
    ##  [13] gdata_2.18.0           magrittr_1.5           checkmate_1.8.5
    ##  [16] memoise_1.1.0          fit.models_0.5-14      cluster_2.0.7-1
    ##  [19] doParallel_1.0.11      fastcluster_1.1.24     annotate_1.58.0
    ##  [22] prettyunits_1.0.2      colorspace_1.3-2       blob_1.1.1
    ##  [25] rrcov_1.4-3            dplyr_0.7.4            RCurl_1.95-4.10
    ##  [28] jsonlite_1.5           genefilter_1.62.0      bindr_0.1.1
    ##  [31] impute_1.54.0          survival_2.42-3        iterators_1.0.9
    ##  [34] glue_1.2.0             hash_2.2.6             gtable_0.2.0
    ##  [37] zlibbioc_1.26.0        XVector_0.20.0         DEoptimR_1.0-8
    ##  [40] scales_0.5.0           mvtnorm_1.0-7          DBI_0.8
    ##  [43] plotrix_3.7            viridisLite_0.3.0      xtable_1.8-2
    ##  [46] progress_1.1.2         htmlTable_1.11.2       foreign_0.8-70
    ##  [49] bit_1.1-12             preprocessCore_1.42.0  Formula_1.2-2
    ##  [52] sqldf_0.4-11           htmlwidgets_1.2        httr_1.3.1
    ##  [55] gplots_3.0.1           RColorBrewer_1.1-2     acepack_1.4.1
    ##  [58] pkgconfig_2.0.1        XML_3.98-1.11          nnet_7.3-12
    ##  [61] locfit_1.5-9.1         dynamicTreeCut_1.63-1  labeling_0.3
    ##  [64] rlang_0.2.0            later_0.7.1            munsell_0.4.3
    ##  [67] tools_3.5.0            gsubfn_0.7             RSQLite_2.1.0
    ##  [70] evaluate_0.10.1        stringr_1.3.0          yaml_2.1.18
    ##  [73] knitr_1.20             bit64_0.9-7            robustbase_0.93-0
    ##  [76] caTools_1.17.1         purrr_0.2.4            bindrcpp_0.2.2
    ##  [79] mime_0.5               compiler_3.5.0         rstudioapi_0.7
    ##  [82] plotly_4.7.1           png_0.1-7              tibble_1.4.2
    ##  [85] geneplotter_1.58.0     pcaPP_1.9-73           stringi_1.1.7
    ##  [88] lattice_0.20-35        Matrix_1.2-14          pillar_1.2.2
    ##  [91] data.table_1.10.4-3    bitops_1.0-6           httpuv_1.4.1
    ##  [94] R6_2.2.2               latticeExtra_0.6-28    promises_1.0.1
    ##  [97] KernSmooth_2.23-15     gridExtra_2.3          codetools_0.2-15
    ## [100] MASS_7.3-50            gtools_3.5.0           assertthat_0.2.0
    ## [103] chron_2.3-52           proto_1.0.0            rprojroot_1.3-2
    ## [106] GenomeInfoDbData_1.1.0 grid_3.5.0             rpart_4.1-13
    ## [109] tidyr_0.8.0            rmarkdown_1.9          shiny_1.0.5
    ## [112] WGCNA_1.63             base64enc_0.1-3

# 4 References

Gu, Z., Liu, J., Cao, K., Zhang, J., & Wang, J. (2012). Centrality-based pathway enrichment: a systematic approach for finding significant pathways dominated by key genes. BMC systems biology, 6(1), 56.

Wu, J. (2016). Transcriptomics and Gene Regulation. Springer.

Zhang, W., Chien, J., Yong, J., & Kuang, R. (2017). Network-based machine learning and graph theory algorithms for precision oncology. npj Precision Oncology, 1(1), 25.