Cross-validation of survival associated biomarkers in gastric cancer using transcriptomic data of 1,065 patients

1Gastric cancer is one of the most common malignancies and displays variable incidence around the globe. About 90–95% of the cases are sporadic, it's incidence is highest in East Asia, Central Eastern Europe, and approximately 75% of cases present in less developed countries. In most developed regions, rates of stomach cancer have variably but uniformly declined over the past decades, fairly due to active surveillance methods in selected societies. There are no solid biomarkers besides HER2 and regular clinicopathological parameters predicting prognosis and response to therapy, and ultimately, there are no efficient therapeutic options available which prove to change the outcome of patients in a groundbreaking manner.

2Following endoscopic examination and histologic confirmation of malignancy in the harvested biopsy, the basis of therapy is still removal of the tumor mass utilizing surgery. The 5-year survival rate for R0 surgical resection ranges from 30 to 50% for patients with stage II and from 10 to 25% for patients with stage III disease. As these patients have a high likelihood of local and systemic relapse, most centers offer them systemic therapy forming the other cornerstone of the treatment. In addition, radiation therapy has proved to improve 5-year survival in resectable tumors.

3Applied prognostic factors of gastric cancer are limited to the clinicopathological properties in the routine setting today, and classically include the WHO histopathological type, Lauren-Järvi classification, size of the tumor, grade, invasion through the gastric wall (pT), vascular invasion, lymph node involvement (pN), etiological background (EBV or Helicobacter pylori) and HER2 overexpression. Biomarkers for diagnosis and prognosis of gastric cancer that have previously been identified are mostly non tumor tissue based, and include carcinoembryonic antigen (CEA), CA 125, CA 19–9, CA 72–4 and alpha-fetoprotein, serum pepsinogen I, and proteases (pepsinogen C, plasminogen activator, matrix metalloproteinases and their inhibitors). Cadherins, mucins and CD44 splicing variants are related to invasion/metastasis and extracellular matrix adhesion and degradation.

4Among tissue based markers, overexpression of human epidermal growth factor receptor 2 (HER2) has been identified as a negative prognostic factor. Trastuzumab with chemotherapy in HER2-positive advanced gastric cancer was investigated in the ToGA study. In this phase 3 trial, 22% of advanced stage cancers overexpressed HER2 and overall survival with trastuzumab was 2.7 months longer (hazard ratio, HR = 0.74, P = 0.0046). In addition, trastuzumab improved all of the secondary end points as well.

5In a search for robust cancer tissue related biomarkers, first we intended to perform a literature review and identify previously described markers for gastric cancer outcome. We merged transcriptomic data of multiple independent datasets to enable a cross-validation of these in a uniform independent cohort. We used uni- and multivariate analyses to assess the prognostic potential for each of the candidate markers. Finally, we compared expression in normal and gastric cancer samples to evaluate the change of the gene expression during tumor formation.

6The entire gastric cancer database includes 1,065 samples from seven independent datasets. Of these, 652 samples were measured with the Affymetrix Human Genome U133 Plus 2.0 Array, 145 with the Human Genome U133A 2.0 Array and 268 with the Human Genome U133A Array. Five arrays did not pass quality control and were excluded from the cross-validation analysis (all five arrays originated in the Bethesda dataset).

7Gender and stage were available for most patients −70% of samples were male and stage III was most common (Figure 1A). Additional clinical parameters including TNM stages, histology and systemic treatment were available for about half the patients – the aggregate clinical characteristics are summarized in Table 1. The median time to first progression (FP) was 18.3 months and the median overall survival (OS) was 28.9 months. Even with these numerically significant differences, the survival curves comparing FP and OS display minor difference (Figure 1B) indicating a short post-progression survival – in the 503 patients with a first event and a known OS, the median post progression survival was 9.4 months.

8Of the clinical parameters, gender, differentiation and histology were not significantly correlated to overall survival. Stage (p = 5.5E-28, see Figure 1C), T (p = 7.9E–15), and N (p = 1.1E-19) delivered high significance while there were not sufficient events to compute correlation to OS for M. Similar results were delivered for FP survival (stage: p = 1.7E-31, T: p = 9.2E-14, and N: p = 4.3E-20). In addition, M was also significant for FP (p = 1.3E-16).

9The keyword search in PubMed resulted in 775 hits, of which 749 were in English language, and 398 were published between 2012–2015. Of these, 40 publications were categorized as review. Following careful and critical evaluation, a list of 29 markers emerged (Supplementary Table 1). Of these candidates, one gene was not present on the gene chips (AFAP1L2), and the remaining 28 were evaluated in the cross-validation.

10Out of the 28 biomarkers, 19 reached significance level with a FDR below 5% for FP and 20 for OS in the univariate analysis investigating gene expression only. Eighteen markers were significant for both FP and OS. Higher expression of BECN1, CASP3, COX2, CTGF, CTNNB1, MET, and SIRT1 correlated to better survival. Higher expression of BIRC5, CNTN1, EGFR, ERCC1, HER2, MMP2, PFKB4, SPHK, SP1, SPARC, TIMP1 and VEGF were predictive for poor outcome. For OS, the direction of correlation to survival was the same for all significant genes. The significant genes with hazard rates and p values are listed in Table 2. Supplementary Table 1. lists the results for all genes.

11Kaplan-Meier survival plots for three of the best performing genes with higher expression correlated to better survival (CASP3, SIRT1 and CTNNB1) and for three of the strongest genes with higher expression correlating to worse survival (BIRC5, HER2 and TIMP- 1) are presented in Figure 2.

12When running a multivariate analysis for OS using all genes, BIRC5 (p = 0.0018), CASP3 (p = 2.3E-04), and CTNNB1 (p = 0.0011) were significant at p < 0.05. Marginal significance was achieved by SP1 (p = 0.07) and VEGF (p = 0.07). When also including stage in the multivariate model, only BIRC5 (p = 0.05) and stage (p = 1.1E-06) were significant. When including stage and age in the multivariate model, only HIF1alpha (p = 0.02), SPARC (p = 0.03), stage (p = 6.8E-07), and age (p = 0.002) were significant.

13In case of FP, when including all genes in a multivariate Cox regression, BIRC5 (p = 0.0017), CASP3 (p = 9.7E-05), CTNNB1 (p = 0.01), MMP-2 (p = 0.0092), SIRT1 (p = 0.035), SPARC (p = 0.0024), and VEGF (p = 0.027) were significant at p < 0.05. However, when including stage or stage and age, only VEGF (p = 0.02), stage (p = 9.3E-07), and age (p = 0.01) remained significant. We have to note that the multivariate analysis used only a fraction of patients included in the univariate analysis, as not all patients had complete clinical annotation (n = 316 for OS and n = 240 for FP).

14We used the expression of MKI67 as a surrogate of proliferation and run a Spearman correlation analysis for all genes. MKI67 itself had a strong prognostic value when examined in a univariate analysis for both overall survival (p = 0.0017, HR = 1.32) and relapse-free survival (p = 0.0015 and HR = 1.39).

15Positive correlation to MKI67 expression was delivered by BIRC5 (coeff = 0.57, p < 1E-20), uPAR (coeff = 0.27, p = 3.9E-19), mTOR (coeff = 0.26, p = 1.6E- 18), SPHK1 (coeff = 0.21, p = 2.1E-12), and HER2 (coeff = 0.21, p = 2.8E-12). Negative correlation was observed for CTGF (coeff = −0.34, p = 9.6E- 31), SPARC (coeff = −0.31, p = 3E-25), PECAM-1 (coeff = −0.30, p = 9.2E-24), and SIRT1 (coeff = −0.23, p = 1.6E-14). As higher expression of multiple genes with negative correlation resulted in better survival (e.g. CTGF, SIRT1), and higher expression of genes with positive correlation delivered worse survival (e.g. BIRC5, HER2, SP1), we computed the correlation between the achieved hazard rate and the correlation coefficient against MKI67 expression. This analysis delivered a borderline significance (coeff = 0.32, p = 0.04). The same analysis performed for HER2 identified SP1 (coeff = 0.26, p = 5.1E-18), BIRC5 (coeff = 0.26, p = 5.2E-18), and EGFR (coeff = 0.20, p = 4.9E-17) having the highest correlation between gene expression and HER2 expression.

16The keywords “gastric” and “normal” GEO delivered 266 datasets. When reducing the search to individual platforms, nine datasets were generated with the GPL96, 35 datasets with the GPL570 and two datasets with the GPL571 platform. Of these, five datasets (GSE44740, GSE51725, GSE13911, GSE43346, and GSE3526) contained expression data for a total of 57 normal gastric tissue specimens.

17When comparing gastric normal and tumor samples, of all 28 genes, 6 were significant below p < 0.01 and had a fold change increase over 1.5 (BIRC5, CTNNB1, HER2, MET, PECAM-1 and uPAR) while only one gene had a 1.5-fold change reduction at the same significance (MMP- 2). The means with 95% confidence intervals for these genes are presented in Figure 3. Supplementary Table 2 contains all the expression values with the Mann-Whitney p value for each gene.

18In this project, we performed a validation for previously identified biomarkers for gastric cancer outcome and prognosis prediction. To assemble a sizeable patient number, we collected data from three cancer centers around the globe (Europe, USA, Australia), and integrated with additional publicly available datasets. Then, in a literature search we selected 28 relevant markers from the past few years (2012–2015), which were disclosed to be related to survival of gastric adenocarcinoma. Analysis was restricted to the most relevant genes by assessing only markers which had been previously published in review manuscripts.

19Most of the molecules were related to worse outcome, being negative prognostic factors. The significant genes include members of the epidermal growth factor family and related genes (EGFR, HER2, SP1), hypoxia-induced and angiogenic genes (HIF1A, VEGF, COX-2), members of the MET and NOTCH signaling pathway (MET, CTNNB1, CNTN1), regulators of survival and proliferation (SP1, CASP3, CTGF), and genes involved in cellular motility (TIMP, MMP2). All these genes are members of important pathways that contribute to progression of the neoplastic process through proliferation and survival by reprogramming the cell's metabolism, mitochondrial functions, protein and lipid synthesis, cytoskeletal organization and signaling.

20The strongest candidate gene was BIRC5 (baculoviral inhibitor of apoptosis repeat containing 5, also known as survivin) – although numerically other genes (like TIMP-1) reached a higher hazard rate, but survivin remained significant in a multivariate analysis as well. BIRC5 has multiple effects including inhibition of apoptosis, enhancing cell proliferation, and promotion of angiogenesis. The correlation between survivin and gastric cancer survival was described in multiple studies. In addition to effect on survival, we observed the highest correlation between BIRC5 and MKI67 expression which supports the link between survivin expression and progression through the cell cycle. However, MKI67 expression per se was not correlated to survival in our Cox regression analysis suggesting that the main effect of BIRC5 on survival is independent of cell proliferation. Theoretically, targeted therapy could be pursued in gastric cancer with siRNA, e.g. LY2181308, is investigated against survivin in multiple types of cancer, however outcome data is pending. Furthermore, immunotherapy may hold promise for these patients as survivin is a strong antigenic peptide for the T cells.

21Among the most promising markers is SIRT1 (Silent mating type information regulation 1), a key gene in the progression of gastric cancer. Earlier, it's expression was associated with histological type, stage, lymph node status, and p53 expression, and proliferation as measured by Ki-67 index. It was also identified as being down-regulated in gastric cancer and a key regulator of NFκB/Cyclin D1 signaling and G1 phase arrest, offering a possible therapeutic intervention in biological models. Here, our results confirm both the association between high SIRT1 and better survival and the correlation between SIRT1 and MKI67 expression. In theory, targeting SIRT1 can be utilized in cancer therapy, mostly cell cycle arrest in G1 phase showed promise in prostate cancer with DU145.

22Previously, receptor tyrosine kinase genomic alterations were detected in 20.6% of cases, affecting ERBB2, FGFR2, and MET, suggesting potential benefit from targeted therapy including MET-amplified gastric tumors and ERBB2 base substitutions. Temporary but durable response to anti-MET agents have already been described. Here, we observed a significant correlation between HER2 and MET and survival, but FGFR2 was not significant. Nevertheless, currently only one of the investigated genes, HER2 expression and amplification is utilized in the routine for prediction of response to anti-HER2 therapy.

23Interestingly, the targetable genes (by administering e.g. cobimetinib, trastuzumab, and ponatinib, respectively, against) MET, HER2, and FGFR2 also displayed the highest difference with a fold change of the mean expressions over 2 when comparing gastric normal and cancer samples. However, again, FGFR2 was not significant – the reason for this is the uneven distribution of FGFR2 expression resulting in almost overlapping median expression at the same time.

24Although our database represents a wide range of clinical cases, the patients are still not sufficiently characterized and this leads to a limitation of our study. While HER2 is an important marker also utilized clinically, our patient samples were collected before the introduction of anti-HER2 therapy. Thus, it was not possible to evaluate the effect of anti-HER2 therapy in the dataset. In the recent years, The Cancer Genome Atlas (TCGA) project proposed a molecular classification dividing gastric cancer into four subtypes. We were also not able to validate the markers in these subtypes as the classification was not available for the patients included in the investigated datasets.

25An alternative approach for survival prediction would be avoiding the utilization of a pre-defined gene to assign patients into prognostic cohorts. Rather, a whole transcriptome gene expression signature could be utilized to select molecularly similar patients and then determine prognostic expectations by evaluating the clinical outcome for these similar patients as has been demonstrated recently. However, no similar methodology has been proposed for gastric cancer, thus we have not included such a model in our meta-analysis.

26In summary, we collected gene expression data sets from three institutions and merged these with public datasets. Then, we performed a literature review and validated previously described markers for gastric cancer outcome. The major advantage of our analysis is that we evaluated all genes in the same set of patients thereby making direct comparison of the markers feasible. The best performing genes include BIRC5, CASP3, CTNNB1, TIMP-1, MMP-2, SIRT, SPARC, and VEGF. The importance of pathological parameters is supported by the fact that only a few genes remained significant when also including stage and age in a multivariate analysis.

27We performed a literature search in PubMed (http://www.ncbi.nlm.nih.gov/pubmed) restricted to the timeframe of 2012 to 2015 utilizing the keywords “gastric”, “cancer”, “survival”, “gene expression” and “biomarker”. To limit the analysis to the most promising markers we selected only the English language manuscripts that were categorized as reviews. We performed the search in November/2015, and then manually continued with review of the publications one by one. Helicobacter pylori infection as a predisposing factor was not analyzed as a prognostic factor of outcome as there was no gene expression marker directly correlated to the infection. We assigned unique gene identifiers for each gene using the HUGO Gene Nomenclature Database (http://www.genenames.org/).

28We assembled a gastric cancer database using samples measured in three different sources including previously partly published data at the Max Delbrück Center for Molecular Medicine, Berlin, Germany (“Berlin dataset”, published in GEO as GSE22377); at the Transgenic Oncogenesis and Genomics Section, Laboratory of Cancer Biology and Genetics, National Cancer Institute, Bethesda, Maryland, USA (“Bethesda dataset”, published in GEO as GSE14210); and at the Peter MacCallum Cancer Centre, Melbourne, Australia (“Melbourne dataset”, published in GEO as GSE51105). Sample collection, hybridization, and gene expression measurements were described previously. The clinical data was updated for each dataset at the end of 2014 and we utilized in the analysis the aggregate database containing all samples with available follow-up data.

29We further extended the database using gene expression data downloaded from GEO. For this, we utilized the keywords “gastric”, “cancer”, “GPL96”, and “GPL570” to search GEO (http://www.ncbi.nlm.nih.gov/geo/). Only publications with available raw data, clinical survival information, and at least 15 patients were included. Affymetrix HG-U133A (GPL96) and HG-U133 Plus 2.0 (GPL570) microarrays were considered because of their overlapping set of 22,277 probe sets and because of our datasets were also derived using these gene chips.

30To discriminate genes related to carcinogenesis, we assembled a database of normal tissues. For this, we used the keywords “gastric” and “normal” in GEO without any limitation regarding publication time or sample number within the study. We included only the GPL96, GPL570, and GPL571 platforms in the search. Samples with premalignant conditions such as intestinal metaplasia were not included as “normal”.

31The raw CEL files were MAS5 normalized in the R statistical environment (http://www.r-project.org) using the Affy Bioconductor library. Quality control for gene chips and control for duplicate samples were performed as described previously. Only arrays passing the quality criteria were utilized. After normalization, only probes measured on both GPL96 and GPL570 were retained (n = 22,277). We subsequently performed a second scaling normalization to set the average expression on each chip to 1000 to reduce batch effects. Kaplan–Meier survival plot and the hazard ratio with 95% confidence intervals and log-rank P values were calculated and plotted in R using Bioconductor packages. False discovery rate (FDR) was computed to correct for multiple testing using the brainwaver library in R as described previously – the FDR cutoff was set at 5%. Expression in cancerous and normal samples was compared using a Mann-Whitney U-test.

32We performed a multivariate analysis using Cox proportional hazards regression including the gene expression markers and clinical variables including stage, age, Lauren classification, differentiation, and gender. In addition to the clinical data, we also determined the HER2 and MKI67 expression using data provided on the gene chips. We computed HER2 status by using the probe set 216836_s_at and setting the cutoff for positivity at 4800. To assess correlation to proliferation, Spearman correlation to MKI67 expression (probes set 212021_s_at) was computed for each of the genes separately. In addition, Spearman correlation was also run for HER2 without using the dichotomization.