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.
Text can be seen as sequential data, either as a sequence of characters or as a sequence of words. Recurrent Neural Network (RNN) operating on sequential data is a type of the neural network. RNN has been applied to a variety of tasks including natural language processing such as machine translation. However, RNN suffers from the problem of long-term dependencies which means that RNN struggles to remember information for long periods of time. An Long Short-Term Memory (LSTM) network is a special kind of RNN that is designed to solve the long-term dependency problem. The bidirectional LSTM network consists of two distinct LSTM networks, termed the forward LSTM and the backward LSTM, which process the sequences in opposite directions. Gated Recurrent Unit (GRU) is a simplified version of LSTM with less number of parameters, and thus, the total number of parameters can be greatly reduced for a large neural network. LSTM and GRU are known to be successful remedies to the long-term dependency problem. The above models take terms of gene sets as input and enrichment scores as output to predict enrichment scores of new terms.
Consider a simple example. Once GSEA is performed, the result calculated from GSEA is fed into the algorithm to train the deep learning models.
library(ttgsea)
library(fgsea)
data(examplePathways)
data(exampleRanks)
names(examplePathways) <- gsub("_", " ", substr(names(examplePathways), 9, 1000))
set.seed(1)
fgseaRes <- fgseaSimple(examplePathways, exampleRanks, nperm = 10000)
data.table::data.table(fgseaRes[order(fgseaRes$NES, decreasing = TRUE),])## pathway pval
## 1: Mitotic Prometaphase 0.0001520219
## 2: Resolution of Sister Chromatid Cohesion 0.0001537043
## 3: Cell Cycle, Mitotic 0.0001255020
## 4: RHO GTPases Activate Formins 0.0001534684
## 5: Cell Cycle 0.0001227446
## ---
## 1416: Downregulation of SMAD2 3:SMAD4 transcriptional activity 0.0022655188
## 1417: HATs acetylate histones 0.0002779322
## 1418: TRAF3-dependent IRF activation pathway 0.0010962508
## 1419: Nephrin interactions 0.0013498313
## 1420: Interleukin-6 signaling 0.0004174494
## padj ES NES nMoreExtreme size
## 1: 0.004481064 0.7253270 2.963541 0 82
## 2: 0.004481064 0.7347987 2.954314 0 74
## 3: 0.004481064 0.5594755 2.751403 0 317
## 4: 0.004481064 0.6705979 2.717798 0 78
## 5: 0.004481064 0.5388497 2.688064 0 369
## ---
## 1416: 0.028982313 -0.6457899 -1.984552 9 16
## 1417: 0.006365544 -0.4535612 -1.994238 0 68
## 1418: 0.017020529 -0.7176839 -2.022102 4 12
## 1419: 0.019558780 -0.6880106 -2.025979 5 14
## 1420: 0.008590987 -0.8311374 -2.079276 1 8
## leadingEdge
## 1: 66336,66977,12442,107995,66442,52276,...
## 2: 66336,66977,12442,107995,66442,52276,...
## 3: 66336,66977,12442,107995,66442,12571,...
## 4: 66336,66977,107995,66442,52276,67629,...
## 5: 66336,66977,12442,107995,66442,19361,...
## ---
## 1416: 66313,20482,20481,17127,17128,83814,...
## 1417: 74026,319190,244349,75560,13831,246103,...
## 1418: 56489,12914,54131,54123,56480,217069,...
## 1419: 109711,14360,20742,17973,18708,12488,...
## 1420: 16194,16195,16451,12402,16452,20848### convert from gene set defined by BiocSet::BiocSet to list
#library(BiocSet)
#genesets <- BiocSet(examplePathways)
#gsc_list <- as(genesets, "list")
# convert from gene set defined by GSEABase::GeneSetCollection to list
#library(GSEABase)
#genesets <- BiocSet(examplePathways)
#gsc <- as(genesets, "GeneSetCollection")
#gsc_list <- list()
#for (i in 1:length(gsc)) {
# gsc_list[[setName(gsc[[i]])]] <- geneIds(gsc[[i]])
#}
#set.seed(1)
#fgseaRes <- fgseaSimple(gsc_list, exampleRanks, nperm = 10000)Since deep learning architectures are incapable of processing characters or words in their raw form, the text needs to be converted to numbers as inputs. Word embeddings are the texts converted into numbers. For tokenization, unigram and bigram sequences are used as default. An integer is assigned to each token, and then each term is converted to a sequence of integers. The sequences that are longer than the given maximum length are truncated, whereas shorter sequences are padded with zeros. Keras is a higher-level library built on top of TensorFlow. It is available in R through the keras package. The input to the Keras embedding are integers. These integers are of the tokens. This representation is passed to the embedding layer. The embedding layer acts as the first hidden layer of the neural network.
if (keras::is_keras_available() & reticulate::py_available()) {
# model parameters
num_tokens <- 1000
length_seq <- 30
batch_size <- 32
embedding_dim <- 50
num_units <- 32
epochs <- 10
# algorithm
ttgseaRes <- fit_model(fgseaRes, "pathway", "NES",
model = bi_lstm(num_tokens, embedding_dim,
length_seq, num_units),
num_tokens = num_tokens,
length_seq = length_seq,
epochs = epochs,
batch_size = batch_size,
use_generator = FALSE,
callbacks = keras::callback_early_stopping(
monitor = "loss",
patience = 10,
restore_best_weights = TRUE))
# prediction for every token
ttgseaRes$token_pred
ttgseaRes$token_gsea[["TGF beta"]][,1:5]
}## pathway
## 614 TGF-beta receptor signaling activates SMADs
## 615 Signaling by TGF-beta Receptor Complex
## 624 Downregulation of TGF-beta receptor signaling
## 1267 TGF-beta receptor signaling in EMT epithelial to mesenchymal transition
## pval padj ES NES
## 614 0.042648445 0.21175102 -0.4407843 -1.551889
## 615 0.004676539 0.04614047 -0.4244875 -1.763638
## 624 0.033777574 0.18377071 -0.4868887 -1.591561
## 1267 0.726840855 0.88913295 -0.2872284 -0.788917Deep learning models predict only enrichment scores. The p-values of the scores are not provided by the model. So, the Monte Carlo p-value method is used within the algorithm. Computing the p-value for a statistical test can be easily accomplished via Monte Carlo. The ordinary Monte Carlo is a simulation technique for approximating the expectation of a function for a general random variable, when the exact expectation cannot be found analytically. The Monte Carlo p-value method simply simulates a lot of datasets under the null, computes a statistic for each generated dataset, and then computes the percentile rank of observed value among these sets of simulated values. The number of tokens used for each simulation is the same to the length of the sequence of the corresponding term. If a new text does not have any tokens, its p-value is not available.
if (exists("ttgseaRes")) {
# prediction with MC p-value
set.seed(1)
new_text <- c("Cell Cycle DNA Replication",
"Cell Cycle",
"DNA Replication",
"Cycle Cell",
"Replication DNA",
"TGF-beta receptor")
print(predict_model(ttgseaRes, new_text))
print(predict_model(ttgseaRes, "data science"))
}## new_text test_value MC_p_value adj_p_value
## 1 Cell Cycle DNA Replication 3.6906314 0.000 0.0000
## 2 Cell Cycle 2.1422198 0.004 0.0080
## 3 DNA Replication 2.5850067 0.002 0.0060
## 4 Cycle Cell 0.7117395 0.210 0.2100
## 5 Replication DNA 1.6805321 0.010 0.0150
## 6 TGF - beta receptor -1.4068508 0.054 0.0648
## new_text test_value MC_p_value adj_p_value
## 1 datum science 0.06976516 NA NAYou are allowed to create a visualization of your model architecture.
## Model: "model"
## ________________________________________________________________________________
## Layer (type) Output Shape Param #
## ================================================================================
## input_1 (InputLayer) [(None, 30)] 0
## embedding (Embedding) (None, 30, 50) 50050
## bidirectional (Bidirectional) (None, 64) 21248
## dense (Dense) (None, 1) 65
## ================================================================================
## Total params: 71,363
## Trainable params: 71,363
## Non-trainable params: 0
## ________________________________________________________________________________Take another exmaple. A set of names of ranked genes can be seen as sequential data. In the result of GSEA, names of leading edge genes for each gene set are given. The leading edge subset contains genes which contribute most to the enrichment score. Thus the scores of one or more genes of the leading edge subset can be predicted.
if (keras::is_keras_available() & reticulate::py_available()) {
# leading edge
LE <- unlist(lapply(fgseaRes$leadingEdge, function(x) gsub(",", "", toString(x))))
fgseaRes <- cbind(fgseaRes, LE)
# model parameters
num_tokens <- 1000
length_seq <- 30
batch_size <- 32
embedding_dim <- 50
num_units <- 32
epochs <- 10
# algorithm
ttgseaRes <- fit_model(fgseaRes, "LE", "NES",
model = bi_lstm(num_tokens, embedding_dim,
length_seq, num_units),
num_tokens = num_tokens,
length_seq = length_seq,
epochs = epochs,
batch_size = batch_size,
verbose = 0,
callbacks = callback_early_stopping(
monitor = "loss",
patience = 5,
restore_best_weights = TRUE))
# prediction for every token
ttgseaRes$token_pred
# prediction with MC p-value
set.seed(1)
new_text <- c("107995 56150", "16952")
predict_model(ttgseaRes, new_text)
}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.
if (keras::is_keras_available() & reticulate::py_available()) {
## 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)
## 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)
## tokenizing text of GSEA
# model parameters
num_tokens <- 5000
length_seq <- 30
batch_size <- 64
embedding_dim <- 128
num_units <- 32
epochs <- 20
# algorithm
ttgseaRes <- fit_model(fgseaRes, "pathway", "NES",
model = bi_lstm(num_tokens, embedding_dim,
length_seq, num_units),
num_tokens = num_tokens,
length_seq = length_seq,
epochs = epochs,
batch_size = batch_size,
callbacks = keras::callback_early_stopping(
monitor = "loss",
patience = 5,
restore_best_weights = TRUE))
# prediction
ttgseaRes$token_pred
set.seed(1)
predict_model(ttgseaRes, c("translation response",
"cytokine activity",
"rhodopsin mediate",
"granzyme",
"histone deacetylation",
"helper T cell",
"Wnt"))
}## R version 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] fgsea_1.26.0 ttgsea_1.8.0 keras_2.11.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 dplyr_1.1.2 tensorflow_2.11.0
## [4] fastmap_1.1.1 textshape_1.7.3 digest_0.6.31
## [7] lifecycle_1.0.3 ellipsis_0.3.2 koRpus_0.13-8
## [10] tokenizers_0.3.0 NLP_0.2-1 magrittr_2.0.3
## [13] compiler_4.3.0 rlang_1.1.0 sass_0.4.5
## [16] tools_4.3.0 utf8_1.2.3 yaml_2.3.7
## [19] qdapRegex_0.7.5 data.table_1.14.8 knitr_1.42
## [22] stopwords_2.3 htmlwidgets_1.6.2 reticulate_1.28
## [25] xml2_1.3.3 textclean_0.9.3 RColorBrewer_1.1-3
## [28] BiocParallel_1.34.0 purrr_1.0.1 grid_4.3.0
## [31] fansi_1.0.4 tm_0.7-11 colorspace_2.1-0
## [34] ggplot2_3.4.2 scales_1.2.1 zeallot_0.1.0
## [37] cli_3.6.1 rmarkdown_2.21 DiagrammeR_1.0.9
## [40] crayon_1.5.2 generics_0.1.3 rstudioapi_0.14
## [43] tfruns_1.5.1 visNetwork_2.1.2 cachem_1.0.7
## [46] sylly.en_0.1-3 parallel_4.3.0 textstem_0.1.4
## [49] base64enc_0.1-3 vctrs_0.6.2 Matrix_1.5-4
## [52] jsonlite_1.8.4 slam_0.1-50 koRpus.lang.en_0.1-4
## [55] lgr_0.4.4 jquerylib_0.1.4 glue_1.6.2
## [58] codetools_0.2-19 cowplot_1.1.1 sylly_0.1-6
## [61] stringi_1.7.12 gtable_0.3.3 munsell_0.5.0
## [64] mlapi_0.1.1 tibble_3.2.1 pillar_1.9.0
## [67] htmltools_0.5.5 float_0.3-1 rsparse_0.5.1
## [70] R6_2.5.1 evaluate_0.20 lattice_0.21-8
## [73] lexicon_1.2.1 png_0.1-8 SnowballC_0.7.0
## [76] syuzhet_1.0.6 RhpcBLASctl_0.23-42 bslib_0.4.2
## [79] text2vec_0.6.3 Rcpp_1.0.10 fastmatch_1.1-3
## [82] whisker_0.4.1 xfun_0.39 pkgconfig_2.0.3Alterovitz, G., & Ramoni, M. (Eds.). (2011). Knowledge-Based Bioinformatics: from Analysis to Interpretation. John Wiley & Sons.
Consoli, S., Recupero, D. R., & Petkovic, M. (2019). Data Science for Healthcare: Methodologies and Applications. Springer.
DasGupta, A. (2011). Probability for Statistics and Machine Learning: Fundamentals and Advanced Topics. Springer.
Ghatak, A. (2019). Deep Learning with R. Springer.
Hassanien, A. E., & Elhoseny, M. (2019). Cybersecurity and Secure Information Systems: Challenges and Solutions and Smart Environments. Springer.
Leong, H. S., & Kipling, D. (2009). Text-based over-representation analysis of microarray gene lists with annotation bias. Nucleic acids research, 37(11), e79.
Micheas, A. C. (2018). Theory of Stochastic Objects: Probability, Stochastic Processes and Inference. CRC Press.
Shaalan, K., Hassanien, A. E., & Tolba, F. (Eds.). (2017). Intelligent Natural Language Processing: Trends and Applications (Vol. 740). Springer.