ttgsea

Introduction

Functional enrichment analysis methods such as gene set enrichment analysis (GSEA) have been widely used for analyzing gene expression data. GSEA is a powerful method to infer results of gene expression data at a level of gene sets by calculating enrichment scores for predefined sets of genes. GSEA depends on the availability and accuracy of gene sets. There are overlaps between terms of gene sets or categories because multiple terms may exist for a single biological process, and it can thus lead to redundancy within enriched terms. In other words, the sets of related terms are overlapping. Using deep learning, this pakage is aimed to predict enrichment scores for unique tokens or words from text in names of gene sets to resolve this overlapping set issue. Furthermore, we can coin a new term by combining tokens and find its enrichment score by predicting such a combined tokens.

Installation

Python is required, although this is a R package. Furthermore, tensorflow and kears should be installed in Python and they should be connected to R, before installing this package (https://tensorflow.rstudio.com, https://keras.rstudio.com).

To install the stable version from Bioconductor: if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("ttgsea") To install the development version from GitHub: devtools::install_github("dongminjung/ttgsea")

Tutorial

The "airway" dataset has four cell lines with two conditions, control and treatment with dexamethasone. By using the package "DESeq2", differntially expressed genes between controls and treated samples are identified from the gene expression data. Then the log2FC is used as a score for GSEA. For GSEA, GOBP for human is obtained from the package "org.Hs.eg.db", by using the package "BiocSet". GSEA is performed by the package "fgsea". Since "fgsea" can accept a list, the type of gene set is converted to a list. Finally, the result of GSEA is fitted to a deep learning model, and then enrichment scores of new terms can be predicted.

``` library(ttgsea)

data preparation

library(airway) data(airway)

differentially expressed genes

library(DESeq2) des <- DESeqDataSet(airway, design = ~ dex) des <- DESeq(des) res <- results(des) head(res)

log2FC used for GSEA

statistic <- res$"log2FoldChange" names(statistic) <- rownames(res) statistic <- na.omit(statistic) head(statistic) ```

The log2FC is calculated for every gene with Ensemble Gene ID. ENSG00000000003 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000938 ENSG00000000971 0.38125398 -0.20681260 -0.03792043 0.08816818 1.37822703 -0.42640213

Using the GOBP gene set, GSEA is performed. ```

gene set

library(org.Hs.eg.db) library(BiocSet) go <- go_sets(org.Hs.eg.db, "ENSEMBL", ontology = "BP") go <- as(go, "list")

convert GO id to term name

library(GO.db) names(go) <- Term(GOTERM)[names(go)]

GSEA

library(fgsea) set.seed(1) fgseaRes <- fgsea(go, statistic) head(fgseaRes[,-"leadingEdge"]) ```

The result of GSEA contains names of pathways and enrichment scores. pathway pval padj log2err 1: 'de novo' AMP biosynthetic process 0.59126984 1.0000000 0.07011322 2: 'de novo' CTP biosynthetic process 0.02202493 0.7889598 0.35248786 3: 'de novo' GDP-L-fucose biosynthetic process 0.72344689 1.0000000 0.06077195 4: 'de novo' IMP biosynthetic process 0.68421053 1.0000000 0.06211242 5: 'de novo' NAD biosynthetic process from tryptophan 0.69094488 1.0000000 0.06211242 6: 'de novo' UMP biosynthetic process 0.79518072 1.0000000 0.05641184 ES NES size 1: 0.5849699 0.9525488 3 2: -0.9523997 -1.4140498 2 3: 0.5771656 0.8519374 2 4: 0.4031617 0.8156349 6 5: 0.4328288 0.8303816 5 6: -0.4293014 -0.7200576 3

The main function of the package "ttgsea" is "fit_model" and it is used to build a deep learning model for predicting new terms. This function gives a trained model. ```

tokenizing text of GSEA

model parameters

numtokens <- 5000 lengthseq <- 30 batchsize <- 64 embeddingdim <- 128 num_units <- 32 epochs <- 20

algorithm

ttgseaRes <- fitmodel(fgseaRes, "pathway", "NES", model = bilstm(numtokens, embeddingdim, lengthseq, numunits), numtokens = numtokens, lengthseq = lengthseq, epochs = epochs, batchsize = batchsize, callbacks = keras::callbackearlystopping( monitor = "loss", patience = 5, restorebestweights = TRUE)) ```

We can monitor the training process. Pearson correlation coefficient is used as a performance metric. Pearson correlation is a measure of how close the predicted value is to the true value. If it is close to 1, the model is considered a good fit. If it is close to 0, the model is not good. A value of 0 corresponds to a random prediction. Epoch 1/20 183/183 [==============================] - 58s 315ms/step - loss: 0.9679 - pearson_correlation: 0.2008 Epoch 2/20 183/183 [==============================] - 60s 326ms/step - loss: 0.8863 - pearson_correlation: 0.3570 Epoch 3/20 183/183 [==============================] - 57s 312ms/step - loss: 0.8132 - pearson_correlation: 0.4489 Epoch 4/20 183/183 [==============================] - 60s 328ms/step - loss: 0.7895 - pearson_correlation: 0.4658 Epoch 5/20 183/183 [==============================] - 58s 319ms/step - loss: 0.7487 - pearson_correlation: 0.5088 Epoch 6/20 183/183 [==============================] - 69s 378ms/step - loss: 0.7373 - pearson_correlation: 0.5322 Epoch 7/20 183/183 [==============================] - 59s 325ms/step - loss: 0.6840 - pearson_correlation: 0.5714 Epoch 8/20 183/183 [==============================] - 59s 321ms/step - loss: 0.6468 - pearson_correlation: 0.6005 Epoch 9/20 183/183 [==============================] - 64s 349ms/step - loss: 0.6341 - pearson_correlation: 0.6109 Epoch 10/20 183/183 [==============================] - 59s 324ms/step - loss: 0.6095 - pearson_correlation: 0.6327 Epoch 11/20 183/183 [==============================] - 59s 321ms/step - loss: 0.5603 - pearson_correlation: 0.6663 Epoch 12/20 183/183 [==============================] - 59s 322ms/step - loss: 0.5289 - pearson_correlation: 0.6937 Epoch 13/20 183/183 [==============================] - 59s 322ms/step - loss: 0.4780 - pearson_correlation: 0.7257 Epoch 14/20 183/183 [==============================] - 59s 323ms/step - loss: 0.4492 - pearson_correlation: 0.7435 Epoch 15/20 183/183 [==============================] - 59s 323ms/step - loss: 0.4181 - pearson_correlation: 0.7654 Epoch 16/20 183/183 [==============================] - 60s 326ms/step - loss: 0.3691 - pearson_correlation: 0.7951 Epoch 17/20 183/183 [==============================] - 59s 324ms/step - loss: 0.3556 - pearson_correlation: 0.8060 Epoch 18/20 183/183 [==============================] - 59s 324ms/step - loss: 1.7955 - pearson_correlation: 0.8179 Epoch 19/20 183/183 [==============================] - 59s 323ms/step - loss: 0.2990 - pearson_correlation: 0.8405 Epoch 20/20 183/183 [==============================] - 60s 327ms/step - loss: 0.2720 - pearson_correlation: 0.8548

The enrichment score of every token is predicted. ```

prediction

ttgseaRes$token_pred ```

The helper T cell seems to have a small negative value. ``` token_term pred 1: SA node 2.914987 2: dicarboxylic acid 2.786132 3: heart induction 2.578254 4: 1 beta 2.493106 5: break process 2.484926

4996: T helper -1.991685 4997: response peptide -2.011833 4998: edit -2.135143 4999: helper T -2.139603 5000: helper cell -2.143030 ```

The trained model is used for predicting new terms. set.seed(1) predict_model(ttgseaRes, c("translation response", "cytokine activity", "rhodopsin mediate", "granzyme", "histone deacetylation", "helper T cell", "Wnt"))

In the table, "testvalue" is the predicted enrichment score. We can get the enrichment score for word or phrase. ``` newtext testvalue MCpvalue adjp_value 1 translation response 1.5176259 0.026 0.06066667 2 cytokine activity 2.1312644 0.000 0.00000000 3 rhodopsin mediate -1.2944517 0.078 0.10500000 4 granzyme -1.5796263 0.014 0.04900000 5 histone deacetylation -0.7845942 0.298 0.29800000 6 helper T cell -1.2319267 0.050 0.08750000 7 Wnt -1.0362499 0.090 0.10500000 ```

The function "plotmodel" draws the plot for model architecture. ``` plotmodel(ttgseaRes$model) ```