Introduction

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 20151. 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 decades2.

Increased rates of early GC detection have increased survival rates for GC patients, but treatment outcomes for GC remain low and difficult to predict3. Moreover, GC is a highly heterogeneous disease as reflected by the numerous histological and molecular classifications4.

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 entities5. After some success in targeted therapies for the treatment of several human cancers6,7, 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 market8,9. Over 90% of drugs fail during the early development stage due to safety concerns or a lack of efficacy10.

The increasing availability of large public datasets such as the Gene Expression Omnibus (GEO) at the National Center for Biotechnology Information (NCBI)11, the Cancer Cell Line Encyclopedia (CCLE)12, DrugBank13, and the Library of Integrated Network-Based Cellular Signatures (LINCS)14,15, which together catalog disease-specific and drug-induced gene expression signatures, offers a time-efficient approach to reposition existing drugs for new indications9,16. 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 drugs17.

Most computational drug repositioning approaches are based on a “guilt by association” strategy18, 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 profiles19. Recently, interest in the use of genomics-based drug repositioning to aid and accelerate the drug discovery process has increased9. 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 drugs20,21.

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

Short Overview of Included Studies

The study selection process is outlined in Fig. 1. 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). Detailed information about the downloaded datasets is summarized in Supplementary Table S2. 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 1
figure 1

Flowchart of the selected process to select gene expression datasets for meta-analysis of gastric cancer. GEO, Gene Expression Omnibus; GC, gastric cancer.

Full size image

Tumor Gene Expression Signatures

The workflow for the exploration of compounds using the calculated RGES values is presented in Fig. 2. 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, log2foldchange >1.2), whereas 53 genes showed decreased expression levels in tumors compared to normal tissues (adjusted P < 0.001, log2foldchange > 1.2; Supplementary Table S3).

Figure 2
figure 2

Workflow 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.

Full size image

Similarity in Gene Expressions between Tumor Samples and GC Cell Lines

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). 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

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 IC50s 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). 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). Likewise, higher compound concentrations (≥10 µM) had lower RGES values than lower concentrations (<10 µM, P < 2.2 × 10−16; Fig. 3C). 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 IC50 values (Spearman correlation rho = 0.3, P = 5.61 × 10−3; Fig. 4).

Figure 3
figure 3

Differences 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.

Full size image
Figure 4
figure 4

Correlation between drug activity and reverse gene expression score (RGES) in AGS cancer cell lines. (A) Correlation between RGES and drug activity (IC50) by linear regression and Spearman’s correlation tests.

Full size image

RGES Summarization and Evaluation

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). 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). 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).

Figure 5
figure 5

Correlation between drug activity (IC50) and summarized reverse gene expression score (sRGES) for all cancer cell lines using a linear regression and a Spearman’s correlation tests.

Full size image
Figure 6
figure 6

Correlation 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.

Full size image

Identification of Reversed Genes and Prediction of Compounds

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). Fifteen compounds, including sorafenib, olaparib, ponatinib, tanespimycin, selumetinib, and elesclomol, were all determined to be active compounds against GC (Supplementary Table S6).

Figure 7
figure 7

Genes 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 IC50. Red and green colors indicate up- and down-regulation, respectively, after compound treatment.

Full size image

Discussion

Methods to identify drug candidates that can reverse the expression states of disease-related genes can complement traditional target-oriented approaches in drug discovery9,22,23,24,25. 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 IC50 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 26 and COL4A1 27 were recently shown to be overexpressed in GC. Meanwhile, the overexpression of UBE2C was related to poor prognosis in GC28 and was a potential biomarker of intestinal type GC29. MIF could also be a potential prognostic factor for GC30. 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 therapy31.

Among the active drugs identified by our analysis, the multiple tyrosine kinase inhibitor sorafenib32 and a poly (ADP-ribose) polymerase (PARP) inhibitor, olaparib33, 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

Collection of Gastric Adenocarcinoma Gene Expression Profiles

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)12.

Gene Expression Data Preprocessing

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 log2 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.34 for meta-analyses of gene expression microarray datasets were followed. The R packages MetaQC35 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

MetaDE36,37,38 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 model39 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

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 Institute40. Cell lines described in LINCS, CCLE, and ChEMBL (version 23 1st May 2017, https://www.ebi.ac.uk/chembl/)41 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

Compound response activity data, described as the half-maximal inhibitory concentrations (IC50) in GC cell lines, were retrieved from ChEMBL. As the IC50 values for a given compound could vary for the same cell line across different studies, the median IC50 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/)42. Sensitivity levels were measured in the form of cellular ATP levels as a surrogate for cell number and growth using CellTiter-Glo assays43. 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 (IC50 < 10 μM) and inactive groups (IC50 ≥ 10 μM) based on their activities in cell lines. An IC50 value of 10 μM was chosen as an activity threshold because compounds with IC50 ≥ 10 μM in primary screenings are often not pursued44.

Reverse Gene Expression Score (RGES) Computation and Summarization

The method used to calculate RGES outcomes was adapted from the previously described Connectivity Map method45. 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 compounds46. 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

In cases for which multiple compound activity IC50 data were available for one compound, median IC50 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 (IC50 < 10 μM) than in the inactive group (IC50 ≥ 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 expressions47. 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

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 IC50 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.