Clinical Significance of Four Molecular Subtypes of Gastric Cancer Identified by The Cancer Genome Atlas Project

Gastric cancer is the fourth most common cancer and third leading cause of cancer-related death worldwide, accounting for an estimated annual 723,100 deaths (1–3). Surgical resection with subsequent adjuvant chemotherapy has been established as an effective treatment for patients with early-stage gastric cancer (4–7). However, recurrence occurs in up to 30% to 40% of patients within 5 years (4, 8–10), suggesting that gastric cancer is a clinically heterogeneous disease. The inherent clinical heterogeneity is most likely due to differences in molecular characteristics of cancer cells. In past years, molecular profiling studies of gastric cancers at a whole-genome level have revealed several genetic and epigenetic changes underlying gastric carcinogenesis (11–13). However, these findings have not been translated to clinical practice, largely owing to the lack of a strong genetic–clinical association.

Recently, as part of The Cancer Genome Atlas (TCGA) project, the genome and proteome of gastric cancer have been extensively characterized to uncover molecular subtypes and identify dysregulated pathways and potential therapeutic targets (14).Integrative analysis of multiple genomic and proteomic data from gastric cancer tissues revealed four molecular subtypes: (i) Epstein–Barr virus (EBV) subtype with extreme DNA hypermethylation, (ii) microsatellite instability (MSI) subtype with elevated mutation rates and hypermethylation, (iii) genomically stable (GS) subtype with less distinctive genomic alterations, and (iv) chromosomal instability (CIN) subtype with marked aneuploidy and frequent focal amplification of receptor tyrosine kinases. Although distinct gastric cancer subtypes with striking genomic features would provide a roadmap for the development of personalized treatment strategies, the clinical relevance of these subtypes, such as their implications for prognosis or response to standard treatment, is currently unknown owing to a lack of sufficient follow-up data from patients in the TCGA cohort. Furthermore, because molecular stratification of gastric cancer in the TCGA study was based on highly complicated integrative analysis of multiple genomic and proteomic data sets, it is very unlikely that these findings will be easily translated into clinical practice.

In the current study, we analyzed gene expression data from TCGA project, uncovered gene expression signatures specific to each of the four molecular subtypes, developed prediction models for stratification of patients with gastric cancer by subtype using these signatures, and tested our model in two large independent cohorts. We also found that the subtypes were predictors of survival outcomes and response to standard adjuvant chemotherapy.

Genomic data from the TCGA gastric cancer cohort were downloaded from the TCGA data portal site (http://cancergenome.nih.gov/) and processed as described in previous studies (14–22). By applying integrative analysis of multiple genomic and proteomic data from gastric cancer tissues, including somatic mutations, mRNA expression, miRNA expression, promoter methylation, somatic copy-number alteration, and protein expression data from reverse phase protein arrays, the TCGA classification scheme employed a decision tree whereby gastric tumors were divided into four subtypes (14). Briefly, tumors were first categorized by the presence of EBV features (EBV subtype), then by the presence of high MSI (MSI subtype). The remaining tumors were further grouped by the number of somatic copy-number alterations: genomically stable (GS subtype) or chromosomal instability (CIN subtype). In the TCGA cohort, mRNA expression data were generated by RNA sequencing, and this information was available for 262 tumor tissue samples. Because most tissues in the TCGA cohort were recently collected, follow-up time for patients in the TCGA cohort was very short and incomplete. Thus, TCGA cohort data from 262 patients were used to generate our model (see below) but were not used for our survival analyses.

Demographic information and clinical data, as well as tissue samples, were obtained from 267 patients with gastric cancer who had undergone gastrectomy as primary treatment between 1999 and 2006 at Korea University Guro Hospital, Kosin University College of Medicine, or Yonsei University Severance Hospital, South Korea. All patients underwent a D2 gastrectomy, and all tissues were snap-frozen and stored in −80°C freezer. Informed consent for sample collection had been obtained from all patients. Generation and analysis of genomic data were carried out at The University of Texas MD Anderson Cancer Center (MDACC, Houston, TX), and these data were designated as the MDACC cohort. Our study was approved by the Institutional Review Boards of MD Anderson Cancer Center and each institute that provided tissues. Of the 267 patients in MDACC cohort, 155 had received standard adjuvant chemotherapy [either single-agent 5-fluorouracil (5-FU) or a combination of 5-FU and cisplatin/oxaliplatin, doxorubicin, or paclitaxel]. Patients in stage I and stage IV with metastasis were not included in subset analysis for assessing benefit of adjuvant chemotherapy.

For an external validation cohort, we used tumor specimen data collected from patients with gastric cancer at the Samsung Medical Center (SMC; n = 432), as described in a previous study (23).

All experiments and analyses were done in the Department of Systems Biology at the MDACC. Gene expression data from the 267 patients in MDACC cohort were generated by hybridizing labeled RNAs to HumanHT-12 v3.0 Expression BeadChips (Illumina). Briefly, total RNA was extracted from the fresh-frozen tissues using a mirVana RNA isolation labeling kit (Ambion). We used 500 ng of total RNA for labeling and hybridization according to the manufacturer's protocols. The microarray data were normalized using the quantile normalization method in the Linear Models for Microarray Data package in the R language environment (24). The expression level of each gene was transformed into a log₂ base before further analysis. Primary microarray data from the MDACC cohort are available in the Gene Expression Omnibus (GEO) database of the National Center for Biotechnology Information (NCBI; accession numbers GSE13861 and GSE26942).

Gene expression data from the SMC cohort were generated using HumanRef-8 WG-DASL v3.0 (Illumina) that contained a subset (24,526 gene features) of probes in Human HT-12. Primary microarray data are available in the GEO database of the NCBI (accession number GSE26253).

Gene expression data were generated and analyzed as explained in previous studies (25–29). BRB-ArrayTools (National Institutes of Health) were used for all statistical analyses of gene expression data (30). We first generated a subtype prediction model using data from the TCGA cohort. For selection of subtype-specific gene sets, multiple two-class t tests were performed for all possible combinations of the four subtypes. Gene expression differences were considered statistically significant if the P value was less than 0.001. Only genes with significant differences in expression in all 3 possible comparisons were considered subtype-specific genes, yielding 349 significant genes for the EBV subtype, 455 for the MSI subtype, 1,513 for the GS subtype, and 143 for the CIN subtype. The top 200 significant genes in each subtype and all 143 genes for the CIN subtype were further selected for development of the prediction model.

To develop a subtype prediction model, we adopted a previously developed model using Bayesian compound covariate predictor algorithms (25–29). Briefly, gene expression data for each subtype gene signature (i.e., the 200 significant genes for each subtype, as described above) were used to generate the Bayesian probability of each tissue sample belonging to a particular subtype. We applied 0.4 as the cutoff of Bayesian probability for each predictor. With this cutoff, the sensitivity and specificity of each predictor ranged from 0.8 to 1 in the training set (the TCGA cohort). Receiver operating characteristic (ROC) analysis of this training set indicated the following order of strength for each predictor: EBV > MSI > GS > CIN (Supplementary Fig. S1); therefore, we adopted a TCGA classification scheme employing a decision tree whereby tumors are grouped into the four subtypes. Briefly, new samples in the test cohorts (i.e., the MDACC and SMC cohorts) were assigned to one of the four subtypes according to Bayesian probability scores. When new samples had more than 2 probability scores above the cutoff value, samples were assigned according to the predetermined strength of the predictors. Samples lacking probability scores above the cutoff value were not assigned to any subtype. Same prediction algorithm was applied to gene expression data from gastric cancer cell lines.

We developed an integrated risk assessment model by pooling the probabilities of the four predictors (subtypes). Because EBV and MSI were associated with good prognosis, we used the inverse of the probability for these subtypes to determine risk of recurrence. GS probability was weighted by a factor of 2 to reflect its strong association with poor prognosis. CIN probability was not modified because it was only moderately associated with poor prognosis.

TCGA Risk Score raw (TRSraw) = (1 − EBV probability) + (1 − MSI probability) + (GS probability × 2) + CIN probability. To generate a dynamic range of scores from 0 to 100, we reformulated TRSraw: TRS = e^TRSraw. This generated TRS values ranging from 3.2 to 85.27. Cutoff points were specified to reflect prognostic differences: low risk (<20), intermediate risk (20–30), and high risk of recurrence (>30).

The association of each subtype with overall survival and recurrence-free survival (RFS) in the MDACC cohort was estimated using Kaplan–Meier plots and log-rank tests. Overall survival was defined as the time from surgery to death, and RFS was defined as the time from surgery to the first confirmed recurrence. Data were censored when a patient was alive without recurrence at last follow-up. Multivariate Cox proportional hazards regression analysis was used to evaluate independent prognostic factors associated with RFS and overall survival, including TRS, tumor stage, and pathologic characteristics as covariates. A P value of less than 0.05 was considered statistically significant. To assess the association of each molecular subtype with benefit from adjuvant chemotherapy, we fitted a Cox proportional hazards model to data from patients in MDACC cohort. All statistical analyses were conducted in the R language environment (http://www.r-project.org). Ingenuity pathway analysis (Ingenuity) was used for gene set enrichment analysis and gene network analysis to identify enriched gene sets and upstream regulators in each subtype.

The subtype-specific gene signatures for each of the four subtypes are shown in Fig. 1 and Supplementary Table S1. We next constructed a prediction model using a Bayesian compound covariate predictor algorithm and tested the strength of each signature (29). When ROC analysis was carried out to the Bayesian probability of each sample, the EBV signature showed highest sensitivity and specificity (100% AUC), and the CIN signature showed the lowest sensitivity and specificity (89.6% AUC). The MSI and GS signatures showed sensitivity and specificity that were lower than EBV but higher than CIN (MSI: 97.7% AUC; GS: 95.4% AUC; Supplementary Fig. S1).

Using the scatter plot matrix approach (Fig. 2), we determined the cutoff of Bayesian probability for each predictor to have reasonably high sensitivity and specificity. Under a cutoff of 0.4, the sensitivity and specificity of each predictor ranged from 0.8 to 1 in the training set (Supplementary Table S2). For construction of the prediction model for the four subtypes, we adopted a TCGA classification scheme that employs a decision tree whereby tumors are grouped into the four subtypes according to Bayesian probability scores: first EBV, then MSI, then GS, and lastly CIN (Fig. 3).

After establishing a robust, highly sensitive and specific prediction model from the TCGA cohort data to categorize patients with gastric cancer by subtype, we examined the association of each subtype with prognosis using gene expression data from the MDACC cohort (n = 267; Supplementary Table S3). When patients in the MDACC cohort were classified by subtype using our prediction model, the EBV subtype was associated with the best prognosis, for both RFS (P = 0.006 by the log-rank test) and overall survival (P = 0.004; Fig. 4A). The GS subtype was associated with the worst prognosis. Patients with the MSI and CIN subtypes had a moderate prognosis that was worse than that of patients with the EBV subtype but better than that of patients with the GS subtype (Fig. 4A).

When the prediction model was applied to the SMC cohort (n = 432; Supplementary Table S3), the GS subtype was associated with the worst prognosis and the EBV subtype was associated with the best prognosis (Fig. 4B). Consistent with the MD Anderson cohort, the MSI subtype was associated with moderate prognosis. However, the prognosis of patients with the CIN subtype was poorer in the SMC cohort than in the MD Anderson cohort, suggesting that the CIN subtype might represent a less homogeneous subgroup.

Because favorable prognosis of patients in EBV subtype is most intriguing, we next carried out gene set enrichment analysis of EBV-specific gene expression signature using ingenuity pathway analysis. It revealed that genes involved in energy production and metabolism are most significantly altered in EBV subtype (Supplementary Table S4). More interestingly, vast majority of metabolic genes are downregulated in EBV subtype (Supplementary Fig. S2).

Patients with the MSI subtype were diagnosed at older ages (median age 60) relative to patients with the other subtypes, whereas those with the GS subtype were diagnosed at relatively younger ages (median age, 52; P = 9.4 × 10⁻⁷ by the Student t test; Supplementary Fig. S3). In addition, most patients with the EBV subtype were male (79%) and diagnosed at younger ages (median age, 53; P = 0.01 by Student t test; Supplementary Fig. S3). This observation is highly consistent with reports from TCGA and previous studies (14, 31, 32), suggesting that our prediction model grouped patients with similar genetic and clinical characteristics to those observed in each subtype in the TCGA data. Taken together, these findings support the robustness of our prediction model and indicate an association between the molecular subtypes and clinical outcomes, suggesting that molecular characteristics of gastric cancer reflected in genomic and proteomic patterns may dictate clinical outcomes.

We next carried out gene network analysis to uncover potential upstream regulators of subtype-specific genes which may contribute to clinical and biological characteristics of each subtype. Not surprisingly, vast majority of predicted upstream regulators of EBV subtype genes are cytokines such as IL1B, IL2, IL3, IL21, IL27, and INFG (Supplementary Table S5), suggesting that infiltrated immune cells in EBV tumors are highly active due to viral infection. Interestingly, in contrast to EBV subtype, many of these cytokines are suppressed in CIN tumors, suggesting that high copy-number alteration may play roles in cancer immunity by suppressing activation of immune cells. This observation is in good agreement with recent study showing strong correlation of tumor aneuploidy with immune evasion across 12 different cancer types (33). MSI subtype is characterized by activation of EZH2, subunit of histone methyltransferase Polycomb repressive complex 2 (PRC2) that functions in the regulation of gene silencing (34), suggesting that further silencing of gene expression by EZH2 may play roles in progression of MSI subtype. Interestingly, GS subtype is characterized by activation of mir-21, a best known oncomir gene (35), as evidenced by downregulation of its many targets genes in GS subtype (Supplementary Fig. S4).

Because adjuvant chemotherapy is the standard treatment for gastric cancer and more than half of patients in the MD Anderson cohort had received adjuvant chemotherapy (5–7), we next sought to determine whether specific subtypes were associated with an increased clinical benefit from adjuvant chemotherapy. We carried out a subset analysis of patients in the MD Anderson cohort with American Joint Commission on Cancer (AJCC) stage II, III, or IV disease without distant metastasis (n = 157); patients with advanced gastric cancer have been shown to benefit most from adjuvant chemotherapy (7). Of 157 patients, 116 received adjuvant chemotherapy. Patients with the CIN subtype exhibited the greatest benefit from adjuvant chemotherapy, as evidenced by significantly increased RFS rates (P = 0.03; Fig. 5A). The 3-year RFS rate was 58.7% for those who received chemotherapy, compared with 33.5% for those who did not. The HR for recurrence among those who received adjuvant chemotherapy was 0.39 [95% confidence interval (CI), 0.16–0.94, P = 0.03; Fig. 5B]. However, no benefit from adjuvant chemotherapy was observed among patients with the GS subtype (P = 0.66; Fig. 5A). The HR for recurrence among those who received adjuvant chemotherapy was 0.83 (95% CI, 0.36–1.89, P = 0.65; Fig. 5B). Patients with the MSI subtype showed only moderate benefit from adjuvant chemotherapy (P = 0.18 by log-rank test; HR, 0.55; 95% CI, 0.22–1.3, P = 0.18; Fig. 5). The benefit of adjuvant chemotherapy could not be assessed for patients with the EBV subtype because all patients received chemotherapy.

To further demonstrate association of subtypes with chemoresistance/chemosensitivity, we used gene expression data and IC₅₀ values of 5-FU from 24 gastric cancer cell lines that are available from Genomics of Drug Sensitivity in Cancer project (36). As seen in gene expression data from tumor tissues, similar fraction of gastric cell lines were grouped into four subtypes according to subtype-specific gene expression signatures (Supplementary Fig. S5). In good agreement with clinical observation, GS cell lines have highest IC₅₀, whereas CIN cell lines have very low IC₅₀, suggesting that GS cells are indeed resistant to 5-FU.

We next developed an integrated risk assessment model (TRS) by pooling the probabilities of the four predictors (subtypes). When patients in the MDACC and SMC cohorts were pooled into a single cohort (n = 699) and stratified according to TRS, the 5-year RFS rate was 66.7% for low-risk patients, 52.1% for intermediate-risk patients, and 37.5% for high-risk patients (P = 5.4 × 10⁻⁶ by log-rank test; Fig. 6A).

To evaluate the prognostic value of the TRS in combination with other clinical variables, we next carried out univariate and multivariate Cox proportional hazards regression analyses with combined clinicopathologic variables in the pooled cohort (MD Anderson + SMC). In addition to T stage, N stage, distant metastasis, and AJCC stage, which are well-known prognostic factors, TRS was a significant predictor of RFS in the univariate analysis (Supplementary Table S6). When we included all relevant clinical variables in a multivariate Cox regression analysis, TRS remained a significant prognostic factor (HR, 1.5; 95% CI, 1.2–1.9, P = 0.001).

Because stage II gastric cancer is considered to be heterogeneous (37), we carried out a subset analysis in patients with stage II disease. In both univariate and multivariate analyses, TRS was a significant predictor of RFS (Supplementary Table S7), providing independent predictive value beyond tumor location, T stage, number of nodes examined, and histologic findings.

TRS showed linear correlation with probability of recurrence at 5 years (Fig. 6B). The distribution of TRS values across the four subtypes reflected the prognostic nature of each subtype; patients with the GS subtype had the highest TRS values and those with the EBV subtype had the lowest TRS values (Supplementary Fig. S6). Taken together, our data strongly suggest that TRS revealed underlying biology not captured by the traditional clinical and pathologic features.

The comprehensive integrative analysis of the genome and proteome of gastric cancer tissues from TCGA uncovered four molecularly distinct subtypes (14).Because the clinical value of subtype classification has been demonstrated in other cancers (38, 39), we developed a genomic prediction model for the subtype classification of gastric cancer using a statistically defined set of multiple genes and gene expression data, and we tested our model in two large independent cohorts (total of 699 patients). This allowed us to properly demonstrate the robustness of the model, as well as the prognostic value of the subtypes in relation to conventional clinical variables.

Our analysis showed that all of the subtypes were present in both test cohorts, clearly demonstrating the reproducibility of tumor classification developed by TCGA and the robustness of the prediction models we developed in the current study. Moreover, clinical characteristics of patients in predicted subtypes are in good agreement with previous observations, further supporting the robustness of our prediction model. More importantly, the four subtypes showed different clinical courses in terms of overall survival and RFS. We found that patients with the EBV subtype had a better prognosis than patients with other subtypes in both cohorts, which is also consistent with a previous report showing that positive immunostaining for EBV was associated with a good prognosis (32, 40, 41). Previous study suggested increased immune response due to viral infection might prevent outgrowth of cancer cells (32).In addition to increased immune response, our data suggested that reduced metabolic activity leading to slow growth of cancer cells might contribute to better prognosis.

Previous TCGA data showed that mutation rates are highest in the MSI subtype and genomic copy-number alterations are highest in the CIN subtype (Supplementary Fig. S7). The EBV subtype is characterized by very high promoter methylation and frequent mutations in PIK3CA, and the GS subtype is characterized by low mutation rates, low copy-number alterations, and frequent mutations in CDH1 and RHOA. Taken together with these TCGA data, our findings suggest that the molecular characteristics in each subtype may dictate clinical outcomes of patients with gastric cancer.

The subtype-specific genetic signatures used in our model predicted not only survival outcomes, but also the relative benefit of adjuvant chemotherapy. Our findings suggested that patients with the CIN subtype most benefitted from adjuvant chemotherapy. Adjuvant chemotherapy was associated with improved outcomes in patients with the CIN subtype but not in patients with the GS subtype. These findings suggest that gastric cancer cells of the GS subtype are resistant to chemotherapy, and frequently altered genes in the GS subtype might account for chemoresistance. For example, NUPR1 is an activated transcription regulator in GS subtype (Supplementary Table S5), and recent studies demonstrated that it enhances chemoresistance in multiple cancers (42–44). However, although our data revealed an association between the CIN subtype and benefit from adjuvant chemotherapy, the retrospective design of our study limited the predictive nature of this association. Therefore, additional studies are warranted to fully assess the ability of this genomic prediction model to identify patients most likely to benefit from adjuvant chemotherapy and those who may need additional treatments.

Using our model, we also developed TRS, an intuitive scoring system ranging from 0 to 100 that can predict risk of recurrence after treatment. Our findings strongly suggest that TRS provides an independent value beyond conventional clinical variables. First, TRS was an independent predictor of recurrence in multivariable analyses, both in the pooled cohort and in the subset of patients with stage II disease. Second, TRS had a strong linear relationship with probability of recurrence at 5 years. Third, the TRS values for each of the subtypes were highly correlated with survival outcomes for the subtypes. Our study has some limitation such as retrospective nature in clinical data and sample collection and lack of diverse ethnic groups in our cohort. Because all patients in current study were from Korea, the clinical association of subtypes needs to be further investigated in larger prospective cohort.

In summary, by analyzing genomic data from TCGA and new data generated in the current study, we demonstrated the clinical significance of the four subtypes of gastric cancer and developed prediction models that can reliably stratify patients with gastric cancer into these four subtypes. Our prediction model could be used to identify not only patients with a poor prognosis (GS subtype), but also those who would most benefit from adjuvant chemotherapy (CIN subtype). Further development of the prediction model will be necessary before it can be implemented into routine clinical practice. Simpler, more robust standardized tools such as qRT-PCR for measuring expression of these genes will be required for clinical use. Increasing the applicability of the prediction model in the clinic will also require selecting fewer genes that well represent each subtype. Nevertheless, the validation of our prediction model in two independent patient cohorts and the fact that the model reflects the biological characteristics associated with each subtype indicate that this prediction model could be used to develop rational therapy recommendations. If confirmed in prospective studies, the association between subtype and adjuvant chemotherapy outcomes might improve patient selection for treatment.

No potential conflicts of interest were disclosed.

Conception and design: H.-S. Lee, W. Jeong, J.-S. Lee

Development of methodology: J.-E. Hwang, H.-S. Lee, J.-S. Lee

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.-E. Hwang, H.-J. Jang, H.-S. Lee, S.C. Oh, J.-J. Shim, S.H. Lee, J.-H. Cheong, W.K. Kang, S. Kim, S.H. Noh, J.-S. Lee

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.-E. Hwang, H.-J. Jang, H.-S. Lee, E.H. Kim, J.-H. Cheong, J. Kim, J. Chae, J.A. Ajani, J.-S. Lee

Writing, review, and/or revision of the manuscript: J.-E. Hwang, H.-J. Jang, H.-S. Lee, J.-J. Shim, K.-W. Lee, J.-H. Cheong, J.A. Ajani, J.-S. Lee

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): E.H. Kim, J.-H. Cheong, J.Y. Cho, S. Kim, J.-S. Lee

Study supervision: H.-S. Lee, W.K. Kang, J.-S. Lee

This work was supported by grants from the NIH (CA127672, CA129906, CA138671, and CA150229); 2016 Institutional Research Grant (IRG) and 2016 Sister Institute Network Fund (SINF) grant from The University of Texas MD Anderson Cancer Center; Korea National Research Foundation (2014M3C1A3051981) funded by the Ministry of Science, ICT & Future Planning; Scientific Research Center Program (2014R1A2A2A01003983); and Research Initiative grant from Korean Research Institute of Bioscience and Biotechnology (KRIBB).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.