Introduction {#Sec1}
============

Gastric cancer (GC) is the fifth most common cancer worldwide and the third leading cause of cancer death, with 1.3 million incident cases and 819,000 deaths occurring globally in 2015^[@CR1]^. Although GC rates have declined in most developed countries, the incidence of non-cardia GC among Caucasians aged 25--39 years has increased in the United States over the past two decades^[@CR2]^.

Increased rates of early GC detection have increased survival rates for GC patients, but treatment outcomes for GC remain low and difficult to predict^[@CR3]^. Moreover, GC is a highly heterogeneous disease as reflected by the numerous histological and molecular classifications^[@CR4]^.

The development of new drugs to treat diseases, especially cancer, is dependent on the identification of novel drug targets. In recent years, an increasing number of innovations have promised to improve our understanding of disease biology, provide novel targets, and catalyze a new era in the development of medicines. However, despite impressive advances in technologies, the situation has remained relatively static in terms of new molecular entities^[@CR5]^. After some success in targeted therapies for the treatment of several human cancers^[@CR6],[@CR7]^, research has focused more on new approaches for the identification of novel targets in cancer therapy. Although large numbers of potential targets have been identified by advanced technologies, it has proven difficult to find targets that are causally involved in the disease.

The number of drugs approved by the US Food and Drug Administration has continuously declined because traditional methods of drug development do not support highly efficient drug discovery. Traditional approaches to develop of new drugs are expensive and time-consuming, with an average of 15 years and a price tag of more than \$2 billion necessary to bring a drug to market^[@CR8],[@CR9]^. Over 90% of drugs fail during the early development stage due to safety concerns or a lack of efficacy^[@CR10]^.

The increasing availability of large public datasets such as the Gene Expression Omnibus (GEO) at the National Center for Biotechnology Information (NCBI)^[@CR11]^, the Cancer Cell Line Encyclopedia (CCLE)^[@CR12]^, DrugBank^[@CR13]^, and the Library of Integrated Network-Based Cellular Signatures (LINCS)^[@CR14],[@CR15]^, which together catalog disease-specific and drug-induced gene expression signatures, offers a time-efficient approach to reposition existing drugs for new indications^[@CR9],[@CR16]^. Several computational methods, such as bioinformatics, system biology, machine learning, and network analysis can be used for drug repositioning or repurposing as well as to identify new indications for drugs^[@CR17]^.

Most computational drug repositioning approaches are based on a "guilt by association" strategy^[@CR18]^, wherein agents having similar properties are predicted to have similar effects. Many drug repositioning strategies are based on different data, including similar chemical structures, genetic variations, and gene expression profiles^[@CR19]^. Recently, interest in the use of genomics-based drug repositioning to aid and accelerate the drug discovery process has increased^[@CR9]^. Drug development strategies based on gene expression signatures are advantageous in that they do not require a large amount of a priori knowledge pertaining to particular diseases or drugs^[@CR20],[@CR21]^.

The purpose of this study is to predict drug candidates that can treat GC using a computational method that integrates publicly available gene expression profiles of GC patient tumors and GC cell lines and cellular drug response activity profiles.

Results {#Sec2}
=======

Short Overview of Included Studies {#Sec3}
----------------------------------

The study selection process is outlined in Fig. [1](#Fig1){ref-type="fig"}. Following the search and selection steps, eight studies: GSE2689, GSE29272, GSE30727, GSE33335, GSE51575, GSE63089, GSE63288, and GSE65801, were included in the final analysis. An additional dataset, GSE54129, was excluded due to lower quantitative QC scores after a MetaQC analysis (Supplementary Table [S1](#MOESM1){ref-type="media"}). Detailed information about the downloaded datasets is summarized in Supplementary Table [S2](#MOESM1){ref-type="media"}. Tumor gene expression signatures were analyzed for 719 GC samples by comparing RNA expression data for 410 tumors and 326 adjacent normal tissues from the GEO. The samples originated from 410 patients, of whom 152 (37.1%) were Korean, 236 (57.6%) were Chinese, and 22 (5.4%) were Caucasians. The samples of patients who had no prior therapy were from GSE29272, GSE65801, and GSE63288. The sample information was not available in GSE30727 nor GSE26899, while the sample information was not mentioned in GSE33335 nor GSE51575. All patients received some type of pre-treatment in GSE63089.Figure 1Flowchart of the selected process to select gene expression datasets for meta-analysis of gastric cancer. GEO, Gene Expression Omnibus; GC, gastric cancer.

Tumor Gene Expression Signatures {#Sec4}
--------------------------------

The workflow for the exploration of compounds using the calculated RGES values is presented in Fig. [2](#Fig2){ref-type="fig"}. All probe sets were re-annotated with the most recent NCBI Entrez Gene IDs and then mapped manually to yield 9,113 unique common genes across the different platforms. A fixed-effect model method was utilized by combining the *P* values in the MetaDE package. Among the gene expression signatures, 136 genes showed increased expression levels in tumors compared to normal tissues (adjusted *P* \< 0.001, log~2~foldchange \>1.2), whereas 53 genes showed decreased expression levels in tumors compared to normal tissues (adjusted *P* \< 0.001, log~2~foldchange \> 1.2; Supplementary Table [S3](#MOESM1){ref-type="media"}).Figure 2Workflow to determine the reverse gene expression score (RGES) using disease and drug gene expression profiles. The public database GEO was used to create cancer gene expression signatures; LINCS L1000 was used as the drug signature database; ChEMBL was used as the drug activity database; and CCLE was used to map cell lines among databases. GEO, Gene Expression Omnibus; GC, gastric cancer; LINCS, The Library of Integrated Network-Based Cellular Signatures; ChEMBL, a chemical database of bioactive molecules maintained by the European Bioinformatics Institute of the European Molecular Biology Laboratory; CCLE, Cancer Cell Line Encyclopedia.

Similarity in Gene Expressions between Tumor Samples and GC Cell Lines {#Sec5}
----------------------------------------------------------------------

The degree of similarity in the gene expressions between tumor samples from the GEO and GC cell lines from the CCLE was assessed by a ranked based Spearman's correlation test. Gene expression profiles for 40 GC cell lines were included in the CCLE (Supplementary Table [S4](#MOESM1){ref-type="media"}). The top 5000 genes in these cell lines were ranked according to their interquartile range across all cell lines used. Among them, less than 0.3% of genes had expression levels in tumor samples from the GEO that did not correlate with those in all GC cell lines. These genes included the ECC10 (0.62%), ECC12 (0.34%), and HGC (0.32%) cell lines (adjusted *P* value \< 0.05).

RGES Computation {#Sec6}
----------------

LINCS data for changes in the expressions of 978 landmark genes after treatment of AGS cell lines with 25 compounds used to treat human gastric adenocarcinoma were used for the RGES computations. The median IC~50s~ values for 2025 compounds used to treat GC cancer cell lines listed in the ChEMBL were used for computation. Disease signatures including 189 DEGs after extraction from the set of LINCS landmark genes were also used for the RGES computation. Variations in the RGES outcomes were evaluated under various biological conditions. The RGES showed larger variations across different cell lines relative to those within different replicates of the same cell line when the same concentration and treatment duration for a compound were used (*P* \< 2.2 × 10^−16^; Fig. [3A](#Fig3){ref-type="fig"}). In addition, longer treatment durations (≥24 h) were associated with lower RGES outcomes compared to shorter durations (\<24 h) when a compound was tested on the same cell line at the same concentration (*P* \< 2.2 × 10^−16^; Fig. [3B](#Fig3){ref-type="fig"}). Likewise, higher compound concentrations (≥10 µM) had lower RGES values than lower concentrations (\<10 µM, *P* \< 2.2 × 10^−16^; Fig. [3C](#Fig3){ref-type="fig"}). The RGES values for the compounds were evaluated by examining correlations with their activities in the same cell line. Finally, the RGES outcomes were correlated with the IC~50~ values (Spearman correlation rho = 0.3, *P* = 5.61 × 10^−3^; Fig. [4](#Fig4){ref-type="fig"}).Figure 3Differences in RGES under the various biological conditions. (**A**) Standard deviation (s.d.) of RGES of individual compounds across different cell lines (grey) vs. across replicates within the same cell line (black grey). (**B**) RGES distribution between treatment durations \<24 hr (grey) and ≥24 hr (black grey). (**C**) RGES distribution between drug concentrations \<10 μM (grey) and ≥10 μM (black grey). Treatment duration and compound concentration were categorized based on compound profiles in LINCS. *P* values waere calculated using a Wilcoxon signed-rank test.Figure 4Correlation between drug activity and reverse gene expression score (RGES) in AGS cancer cell lines. (**A**) Correlation between RGES and drug activity (IC~50~) by linear regression and Spearman's correlation tests.

RGES Summarization and Evaluation {#Sec7}
---------------------------------

Summarized RGES (sRGES) values were computed by weighting various cell lines, compound concentrations, and treatment durations. A number of known methods were used to summarize the RGES and obtain sRGES values (Supplementary Table [S5](#MOESM1){ref-type="media"}). The calculated sRGES scores for each compound were significantly correlated with drug activity (Spearman correlation rho = 0.27 and *P* = 1.04 × 10^−2^; Fig. [5](#Fig5){ref-type="fig"}). Additionally, CTRP was used as an external dataset to confirm the correlation between reversal potency and compound activity. Activity data expressed as AUC values for 546 compounds tested in GC cell lines were collected from CTRP. After the sRGES computation, the median AUC values across multiple cell lines were used to evaluate the sRGES. The sRGES values were significantly correlated with the AUC values (rho = 0.368, *P* = 3.8 × 10^−8^; Fig. [6](#Fig6){ref-type="fig"}).Figure 5Correlation between drug activity (IC~50~) and summarized reverse gene expression score (sRGES) for all cancer cell lines using a linear regression and a Spearman's correlation tests.Figure 6Correlation between AUC and sRGES of compounds. AUC data were retrieved from CTRP. Sensitivity were measured in the form of cellular ATP levels as a surrogate for cell number and growth using a CellTiter-Glo assay. A compound-performance score was computed at each concentration of compound. The AUC using percent-viability scores was computed as a metric of sensitivity. Median values were used to summarize AUC across all gastric cancer cell lines examined. A Spearman's correlation test was used to analyze the correlation between sRGES and AUC. AUC, areas under the concentration-response curve.

Identification of Reversed Genes and Prediction of Compounds {#Sec8}
------------------------------------------------------------

Using the correlation between the sRGES outcome and the compound activity, compounds having high reversal potency for GC were identified. Next, genes having expression levels that were reversed by the active compounds were predicted by a leave-one-compound-out approach. The five genes that showed significant reversals of expression following treatment with the GC cell lines with the active compounds included: (i) the collagen type IV α1 chain (*COL4A1*); (ii) procollagen-lysine 2-oxoglutarate 5-dioxygenase 3 (*PLOD3*); (iii) ubiquitin conjugating enzyme E2 C (*UBE2C*); (iv) macrophage migration inhibitory factor (*MIF*); and (v) pre-mRNA processing factor 4 (*PRPF4*) (Fig. [7](#Fig7){ref-type="fig"}). Fifteen compounds, including sorafenib, olaparib, ponatinib, tanespimycin, selumetinib, and elesclomol, were all determined to be active compounds against GC (Supplementary Table [S6](#MOESM1){ref-type="media"}).Figure 7Genes having reversed expression in response to treatment with active and inactive  compounds. Low rank and high rank suggests that the gene expression is down- and upregulated, respectively, by the corresponding compound. The heatmap indicates the relative position of a gene in a ranked drug expression profile. Position are normalized and compound columns are ordered according to IC~50~. Red and green colors indicate up- and down-regulation, respectively, after compound treatment.

Discussion {#Sec9}
==========

Methods to identify drug candidates that can reverse the expression states of disease-related genes can complement traditional target-oriented approaches in drug discovery^[@CR9],[@CR22]--[@CR25]^. In this study, we used public cancer genomic and pharmacologic databases to demonstrate the reversal potency relationship between DEGs and drug activity and to predict potential new drug candidates for GC.

Our results showed that the ability of drugs to reverse DEGs was correlated with drug activity in GC, although this correlation was highly dependent on the cell line as well as the drug concentration and treatment duration. The positive correlation between sRGES and IC~50~ values indicated that combining disease gene expression data derived from clinical samples with drug gene expression profiles obtained from results with *in vitro* cell lines could be used to predict drug activity.

In our study, five GC genes, *COL4A1*, *PLOD3*, *UBE2C*, *MIF*, and *PRPF4*, showed reversed expressions in response to 15 active compounds. To the best of our knowledge, this is the first study of drug repositioning using a computational reversal gene expression approach in GC. Among these genes, *PLOD3*^[@CR26]^ and *COL4A1*^[@CR27]^ were recently shown to be overexpressed in GC. Meanwhile, the overexpression of *UBE2C* was related to poor prognosis in GC^[@CR28]^ and was a potential biomarker of intestinal type GC^[@CR29]^. *MIF* could also be a potential prognostic factor for GC^[@CR30]^. These genes showed reversed expression levels and thus may be feasible as therapeutic targets for GC. Additionally, *PRPF4* as a pre-mRNA splicing factor has been suggested as a potential therapeutic target for cancer therapy^[@CR31]^.

Among the active drugs identified by our analysis, the multiple tyrosine kinase inhibitor sorafenib^[@CR32]^ and a poly (ADP-ribose) polymerase (PARP) inhibitor, olaparib^[@CR33]^, have completed phase II and phase III clinical trials, respectively, for GC patients. Meanwhile, the heat shock protein 70 (Hsp70) inducer elesclomol, the novel tyrosine kinase inhibitor ponatinib, the heat shock protein 70 (Hsp90) inhibitor tanespimycin, and the mitogen-activated protein kinase inhibitor selumetinib have not been previously studied clinically for their effectiveness against GC.

GC is a heterogeneous disease that involves multiple factors associated with various molecular pathways that can function differently during the cancer development process. A limitation of this study is that the GC disease gene expression datasets from the GEO are not uniformly associated with clinical outcomes or GC etiologies. The drug activity of predicted compounds may also vary because the GC disease states varied for individual patients. Sampling time information is important, as samples obtained after the initial neo-adjuvant chemotherapy can affect the results of this meta-analysis. Nonetheless, such information was not available from some datasets.

Many recent projects focus on precision medicine to provide insights between diseases and genes. A repurposing strategy based on alterations of driver genes in each tumor can be used to identify therapeutic targets. The collection of therapeutic agents targeting driver genes and determining the connection between each patient and the targeted therapies can enhance promising drug repositioning opportunities and eventually benefit patients. Therefore, RGES may improve predictions of drug candidates because it is based on the molecular characteristics of actual tumors.

Therapeutic efficacy is more complex than a simple correlation of gene expression profiles with drugs and diseases. Therefore, our findings with regard to drug candidates will require further preclinical testing and demonstrations in clinical trials, although our results did validate that the method of the computational analysis of public gene expression databases is a potentially useful means of drug discovery. In summary, our computational approach combined disease gene expression with drug-induced expression profiles in GC to identify new drugs and target genes for GC therapy. This approach can also be used to predict the efficacy of new drug candidates with which to treat GC. This computational approach could be broadly applied to other diseases for which reliable gene expression data are available.

Materials and Methods {#Sec10}
=====================

Collection of Gastric Adenocarcinoma Gene Expression Profiles {#Sec11}
-------------------------------------------------------------

Publicly available gene expression profiles for GC patients were downloaded from the GEO database of the NCBI (<https://www.ncbi.nlm.nih.gov/geo>), A search of the GEO database was conducted in July of 2018 using 'gastric cancer' as a key search phrase. The results for deposits made since January of 2015 were filtered using the search terms *Homo sapiens*, expression profiling by array, and expression profiling by high-throughput sequencing. Only original experimental datasets that compared the expression levels of mRNAs between GC tumors and normal tissue controls were selected. Datasets containing more than ten sets of normal and tumor samples were retained. Additionally, gene expression profiles of human gastric adenocarcinoma cell lines were downloaded from the CCLE (version 2.7. updated 2015 <https://portals.broadinstitute.org/ccle>)^[@CR12]^.

Gene Expression Data Preprocessing {#Sec12}
----------------------------------

The GEO accession number, platform, sample type, numbers of cases and controls, references and expression data were extracted from each of the identified datasets, which were then individually preprocessed using a log~2~ transformation and normalization approach. If there were multiple probes for the same gene, the probe values were averaged for that gene expression level. All probe sets on different platforms were re-annotated to use the most recent NCBI Entrez Gene Identifiers (Gene IDs), and the Gene IDs were used to cross-map genes among the different platforms. Only genes present in all selected platforms were considered. To combine the results from individual studies and to obtain a list of more robust DEGs between GC and normal control tissues, guidelines outlined by Ramasamy *et al*.^[@CR34]^ for meta-analyses of gene expression microarray datasets were followed. The R packages MetaQC^[@CR35]^ was used for quality control (QC). MetaQC uses six quantitative QC parameters: (i) measures of internal QC; (ii) measures of external QC; (iii) accuracy QC of featured genes; (iv) accuracy QC of the pathway; (v) consistency of QC in the ranking of featured genes; and (vi) consistency QC in the ranking of the pathway. The mean rank of all QC measures in each dataset was also determined as a quantitative summary score by calculating the ranks of each QC measure among all included datasets.

Disease Gene Expression Signatures {#Sec13}
----------------------------------

MetaDE^[@CR36]--[@CR38]^ was used to the identify DEGs in GC. A moderated *t*-statistic was used to calculate the *P* values for each dataset, and a meta-analysis was conducted with a fixed-effect model^[@CR39]^ using the MetaDE package to identify DEGs. Additionally, similarities among gene expressions profiles between tumor samples from the GEO and GC cell lines from the CCLE were assessed.

Compound Gene Expression Profiles {#Sec14}
---------------------------------

Level 4 gene expression profiles consisting of 978 landmark genes (L1000 genes) from LINCS as of May of 2018 were downloaded from LINCS cloud storage (<https://lincscloud.org/>), hosted by the Broad Institute^[@CR40]^. Cell lines described in LINCS, CCLE, and ChEMBL (version 23 1st May 2017, <https://www.ebi.ac.uk/chembl/>)^[@CR41]^ were mapped using GC cell line names followed by manual inspections. Meta-information for compound-induced gene expressions, including the cell line types as well as the treatment durations and drug concentrations was retrieved. Only small-molecule perturbagens having high-quality gene expression profiles (is_gold = 1, annotated in the meta-information) were used for further analysis.

Compound activity profiles {#Sec15}
--------------------------

Compound response activity data, described as the half-maximal inhibitory concentrations (IC~50~) in GC cell lines, were retrieved from ChEMBL. As the IC~50~ values for a given compound could vary for the same cell line across different studies, the median IC~50~ value was used. Compounds included in the ChEMBL and LINCS were mapped using The International Union of Pure and Applied Chemistry International Chemical Identifier keys. Additionally, the area under the curve (AUC) values for compound activity data in GC cell lines were retrieved from the Cancer Therapeutic Response Portal (CTRP ver 2, <https://portals.broadinstitute.org/ctrp.v2.1/>)^[@CR42]^. Sensitivity levels were measured in the form of cellular ATP levels as a surrogate for cell number and growth using CellTiter-Glo assays^[@CR43]^. A compound-performance score was computed at each concentration of compound. The AUC using percent-viability scores was computed as a metric of sensitivity given that AUC reflects both relative potency and the total level of inhibition observed for a compound across CCLs. Median AUC values across multiple cell lines were used. Compounds were categorized into active (IC~50~ \< 10 μM) and inactive groups (IC~50~ ≥ 10 μM) based on their activities in cell lines. An IC~50~ value of 10 μM was chosen as an activity threshold because compounds with IC~50~ ≥ 10 μM in primary screenings are often not pursued^[@CR44]^.

Reverse Gene Expression Score (RGES) Computation and Summarization {#Sec16}
------------------------------------------------------------------

The method used to calculate RGES outcomes was adapted from the previously described Connectivity Map method^[@CR45]^. Briefly, genes were initially ranked by their expression values for each compound. An enrichment score for each set of upregulated and downregulated disease genes was computed based on the positions of the genes in the ranked list. RGES values emphasize the reversal correlation by capturing the reversal relationship between the DEGs and compound-induced changes in gene expressions. Therefore, a lower negative RGES indicates a greater likelihood of reversing changes in disease gene expressions, and *vice versa*. In addition, Spearman's correlation coefficient, Pearson correlation coefficient, and cosine similarity were computed between the DEGs and compound activities as an alternate means of computing the reversal relationship between DEGs and active compounds^[@CR46]^. The databases can list multiple gene expression profiles associated with one compound due to testing in various cell lines, compound treatment concentrations, and compound treatment durations, which resulted in multiple RGES outcomes for one compound that could reverse disease gene expression. Given these variations, sRGES were weighted and calculated. Results obtained for 10 μM drug concentrations and 24 h treatments were used to define the reference conditions. The analysis code and an example are provided at <https://github.com/Bin-Chen-Lab/RGES>.

Identification of Reversed Genes {#Sec17}
--------------------------------

In cases for which multiple compound activity IC~50~ data were available for one compound, median IC~50~ values were calculated. In cases for which multiple gene expression profiles yielded multiple RGES values for one compound, a median RGES value was calculated from the GC cell lines. Each gene expression profile was sorted according to its expression value. Upregulated genes were ranked high (i.e., on the top), whereas downregulated genes were ranked low (i.e. on the bottom). Among the upregulated genes, reversal genes were defined as those that were ranked lower in the inactive group (IC~50~ \< 10 μM) than in the inactive group (IC~50~ ≥ 10 μM). In contrast, among the downregulated genes, the reversal genes were defined as those that were ranked higher in the active group than in the inactive group. A leave-one-compound-out cross-validation approach was used to find genes having reversed expressions^[@CR47]^. For each trial, one compound was removed and the reversed genes were then identified using the approach described above. Only those genes that were significantly reversed in all trials were retained. Genes having *P* \< 0.1 in all trials were considered as reversal genes.

Statistical Analysis {#Sec18}
--------------------

The degrees of similarity in the gene expressions between tumor samples from the GEO and GC cell lines from the CCLE were assessed by Spearman's rank correlation testing, as was the similarity of RGES and IC~50~ from ChEMBL or AUC from CCLE. A Wilcoxon signed-rank test was used to assess differences between RGES across the same and different cell lines, longer (fferent cell lines, \<24 h) treatment durations, and higher (≥10 μM) and lower (\<10 μM) drug concentrations. *P* values were adjusted with a Benjamini and Hochberg's false discovery rate method to correct for multiple testing.

Supplementary information
=========================

 {#Sec19}

Supplementary materials

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

In-Wha Kim and Jung Mi Oh Contributed equally.

Electronic supplementary material
=================================

**Supplementary information** accompanies this paper at 10.1038/s41598-019-39228-9.

This study was supported by the National Research Foundation of Korea grant funded by the Korea government (MSIT) (No. NRF-2018R1A2B6001859).

I.W.K. and J.M.O. created and designed the study. I.W.K., J.H.K. and H.J. contributed to data collection and analyzed the data. I.W.K. contributed to write manuscript. J.M.O. and M.G.K. reviewed the manuscript. S.K. reviewed the method and the manuscript. All authors read and approved the final manuscript. J.M.O. and I.W.K. contributed equally to this work.

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing Interests {#FPar1}
===================

The authors declare no competing interests.
