A five-gene signature for predicting overall survival of esophagus adenocarcinoma.

ABSTRACT: Esophageal adenocarcinoma (EAC) is common and aggressive with increasing trend of incidence. Urgent need for an effective signature to assess EAC prognosis and facilitate tailored treatment is required.Differentially expressed mRNAs (DEMs) were identified by analyzing EAC tissues and adjacent normal samples from The Cancer Genome Atlas (TCGA). Then univariate regression analyses were performed to confirm prognostic DEMs. We used least absolute shrinkage and selection operator (LASSO) to build a prognostic mRNA signature whose performance was assessed by Kaplan-Meier curve, receiver operating characteristic (ROC). GSE72874 were used as an external test set. The performances of the signature were also validated in internal TCGA and external test sets. Gene set enrichment analysis (GSEA) and tumor immunity analysis were performed to decipher the biological mechanisms of the signature.A 5-mRNA signature consisted of SLC26A9, SINHCAF, MICB, KRT19, and MT1X was developed to predict prognosis of EAC. The 5-mRNA signature was promising as a biomarker for predicting 3-year survival rate of EAC in the internal test set, the entire TCGA set, and the external test set with areas under the curve (AUC) = 0.849, 0.924, and 0.747, respectively. Patients were divided into low- and high-risk groups based on risk scores of the signature. The high-risk group was mainly associated with cancer-related pathways and low levels of B cell infiltration.The 5-mRNA prognostic signature we identified can reliably predict prognosis and facilitate individualized treatment decisions for EAC patients.

Introduction

Esophagus cancer (ESCA) is the 6th most lethal cancer and the 8th most prevalent malignancy globally.[ The 2 major subtypes of ESCA include esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC). In recent decades, ESCC has decreased in prevalence whilst EAC has gradually increased.[ EAC is a particularly aggressive form of ESCA with overall 5-year survival rates as low as ∼18%.[ The lack of early stage diagnostics coupled to high rates of metastasis and drug resistance contribute to the poor outcomes of EAC.[ Despite the improvements in EAC therapeutics, their benefits must be balanced with their side effects. Precision therapy based on the molecular basis for EAC development offers the most hope for effective therapeutic interventions. Although tumor–node–metastasis (TNM) staging has been applied for prognostic prediction and individual treatment,[ the patients with the same stage can have significantly different outcomes in clinical practice. As such, urgent need for effective models to assess EAC prognosis and facilitate tailored treatment are required.

In recent years, with development of microarray and high-throughput sequencing, studies have focused on mRNAs for the prediction of EAC prognosis. Kim et al identified a 2-gene signature (SPARC and SPP1) associated with prognosis in EAC.[ In addition, dysregulated DCK, PAPSS2, SIRT2, and TRIM44 have been linked to overall survival. As such, a 4-gene signature based on these genes was independently developed for EAC prognostics.[ By evaluating the gene expression profiles in 64 patients with EAC, Pennathur et al constructed an internally cross-validated 59-mRNA prognostic signature.[ There were however limitations in these studies, including a lack of the cross-platform validation datasets, small sample sizes, and the absence of external testing.

The Cancer Genome Atlas (TCGA) is a cancer genome project that provides large-scale genomic and clinical information spanning 33 cancer types. The Gene Expression Omnibus archives Array- and sequence-based data.[ Here, we explored the available EAC data from the TCGA and GEO, and developed a 5-mRNA signature for EAC. This signature was validated in independent external test set, highlighting its promise for EAC prognostics.

Material and methods

Data sources and processing

RNAseq data and clinical data were downloaded from TCGA database (https://portal.gdc.cancer.gov/). Gene annotations were obtained from GENCODE datasets (www.gencodegenes.org). We converted the number of fragments per kilobase of non-overlapped exons per million fragments mapped into the value of transcripts per kilobase million (TPM) to mimic those of the microarrays.[ Only mRNAs with TPM values >0.1 in at least 50% of the ESCA samples were enrolled for analysis. Expression was defined as log2(x + 1), whilst x represented the TPM value. To predict the prognosis of EAC, we excluded samples with a <30-day censor time. We obtained 86 samples from the TCGA including 75 EAC cases and 11 normal tissue samples. Differentially expressed mRNAs (DEMs) were identified using the R package “limma” in EAC versus normal tissue. Selection criteria for significant DEMs were P < .05 and |log2 fold change| >1.5.

We downloaded GSE72874 from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database as an external validation set using R package “GEOquery.” Corresponding clinical data were achieved from Supplemental data submitted by Krause et al.[ Since the data came from the TCGA database and the Gene Expression Omnibus database, no ethical approval was required.

The tumor immune estimation resource (TIMER) algorithm (https://cistrome.shinyapps.io/timer/) was used to measure the abundance of B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells in EAC.[ This provided details on the immune infiltration in EAC and its relevance to the identified gene panel signature.

Construction and validation of the mRNA signature

To further screen DEM correlating with prognosis, univariate COX proportional hazards regression analysis was performed. Only prognostic DEMs with P < .05 were considered as statistically significant mRNAs.

We divided the TCGA datasets into training (n = 37) and tests (n = 38). Least absolute shrinkage and selection operator (LASSO) regression was selected as the optimal dimension reduction technique,[ previous reported studies have applied it to filter significant genes.[ LASSO was performed to screen and confirm the selected DEMs in the training set. Risk scores for each patient were calculated with the linear combinational of the selected gene expression weighted by their regression coefficients. The formula was set as follow:

Risk score = β1 × gene 1 + β2 × gene 2 + … + βn × gene n,

where β is the coefficient of each gene, and gene indicates relative expression of gene.

Datasets were divided into low- and high-risk groups using the median score. Survival curves between these groups were assessed via the Kaplan–Meier method, whilst P-values were estimated using log-rank tests. Receiver operating characteristic (ROC) and areas under the curve (AUC) were calculated to determine the accuracy of the mRNA signature by using the R “timeROC” package. Further, validations in the internal and external tests were performed using the risk score formula.

Biological functions of the mRNA signature

Gene set enrichment analysis (GSEA) was performed using the R package “clusterProfiler”[ to study the biological mechanisms between low- and high-risk EAC groups. Kyoto Encyclopedia of Genes and Genomes (KEGG, c2.cp.kegg.v7.0.entrez.gmt) and Reactome (c2.cp.reactome.v7.0.entrez.gmt) datasets were selected from the molecular signature database. Gene Ontology (GO)[ of the Biological Process (GO_BP), Cellular Components (GO_CC), and Molecular Functions (GO_MF) was performed to examine the biological functions of the signature. Cut-off values for significant enriched terms were P < .05.

Statistical analyses

The R statistical package (R version 3.5.2; http://www.Rproject.org) was used to perform all statistical analyses and for plot construction. We used “survival” package of R to carry out univariate and multivariate Cox proportional hazard assessments. Hazard ratios (HR), 95% confidence intervals, and P-values were calculated. A Student t test was performed to compare continuous variables. A 1-way analysis of variance was used for multiple group comparisons. The association between tumor immunity infiltration and gene expression level was estimated using the Spearman correlation test. P-values < .05 were considered significant.

Results

Study datasets

A flowchart of the study is shown in Figure 1. We obtained data from 75 patients with EAC from the TCGA dataset, divided into training (n = 37) and internal test sets (n = 38). In total, 43 patients with EAC from the GSE72874 dataset were enrolled for further analysis. Basic characteristics in these 2 datasets are shown in Table 1.

Figure 1

A flow diagram of establishment and validation of the prognostic mRNA signature. DEMs = differentially expressed mRNAs, EAC = esophageal adenocarcinoma, LASSO = least absolute shrinkage and selection operator, TCGA = The Cancer Genome Atlas.

Table 1

Clinical characteristics of patients with EAC in TCGA and GSE72874 datasets.

	TCGA		GSE72874
	Training dataset	Testing dataset	External validation dataset
Sample number	37	38	43
Survival time (yr [SD])	1.29 (1.24%)	1.57 (1.44%)	1.71 (1.36%)
Status
Alive	23 (62.2%)	19 (50.0%)	26 (60.5%)
Dead	14 (37.8%)	19 (50.0%)	17 (39.5%)
Age
<60	10 (27.0%)	18 (47.4%)	11 (25.6%)
>60	27 (73.0%)	20 (52.6%)	32 (74.4%)
Sex
Female	5 (13.5%)	6 (15.8%)	3 (7.0%)
Male	32 (86.5%)	32 (84.2%)	40 (93.0%)
Height (cm)
<175	18 (48.6%)	17 (44.7%)
>175	16 (43.2%)	19 (50.0%)
NA	3 (8.1%)	2 (5.3%)
Weight (kg)
<85	19 (51.4%)	24 (63.2%)	40 (93.0%)
>85	17 (45.9%)	14 (36.8%)	3 (7.0%)
NA	1 (2.7%)	0 (0.0%)	0 (0.0%)
Race
Non-white	0 (0.0%)	1 (2.6%)
White	28 (75.7%)	30 (78.9%)
NA	9 (24.3%)	7 (18.4%)
Alcohol history
No	14 (37.8%)	12 (31.6%)
Yes	22 (59.5%)	26 (68.4%)
NA	1 (2.7%)	0 (0.0%)
Barrett disease
No	21 (56.8%)	25 (65.8%)
Yes	14 (37.8%)	10 (26.3%)
NA	2 (5.4%)	3 (7.9%)
Tumor size
I+II	14 (37.8%)	15 (39.5%)
III+IV	22 (59.5%)	22 (57.9%)
NA	1 (2.7%)	1 (2.6%)
Node status
Negative	11 (29.7%)	8 (21.1%)
Positive	24 (64.9%)	28 (73.7%)
NA	2 (5.4%)	2 (5.3%)
Metastasis
0	29 (78.4%)	26 (68.4%)
1	2 (5.4%)	8 (21.1%)
NA	6 (16.2%)	4 (10.5%)
Stage
I+II	16 (43.2%)	16 (42.1%)
III+IV	21 (56.8%)	20 (52.6%)
NA	0 (0.0%)	2 (5.3%)

EAC = esophageal adenocarcinoma, NA = not available, SD = standard deviation, TCGA = The Cancer Genome Atlas.

DEMs associated with EAC prognosis

Through the analysis of EAC and normal tissues in the TCGA dataset, we extracted 852 upregulated and 160 downregulated genes (Supplemental Fig. S1). By subjecting DEMs to univariable COX regression analysis, 138 differentially expressed genes (DEGs) related to overall survival were identified for the prognostic signature (Supplemental Table S1).

Construction of the mRNA signature and evaluation of its prognostic ability in the training set

Lambda value was set using the lambda.min. Five mRNAs with non-zero coefficients were identified (Supplemental Fig. S2). We created a risk score formula as follows:

Risk score =

(−0.0755 × relative expression of SLC26A9) +

(0.6688 × relative expression of SINHCAF) +

(0.3071 × relative expression of MICB) +

(−0.5899 × relative expression of KRT19) +

(−0.4242 × relative expression of MT1X)

In the training set, we classified patients into high (n = 19) and low risk groups (n = 19) using a median risk score of −1.60 as the cut-off. Details on the survival status according to risk scores are shown in Figure 2A. Upon analysis of the heatmaps of the mRNAs (Fig. 2D), 2 mRNAs had positive coefficients including SINHCAF and MICB which indicated that the higher expression levels of these mRNAs were linked with poor survival. The other 3 mRNAs with negative coefficients were SLC26A9, KRT19, and MT1X which demonstrated that it was a positive correlation between their expression levels and clinical outcome. Survival curves indicated that patients in the lower-risk group had higher survival rates (Fig. 2B). Time-dependent ROC analyses demonstrated that AUCs of the mRNA-based signature were 0.760, 0.807, and 0.849 at 1-, 2-, and 3-year survival times, respectively (Fig. 2C).

Figure 2

Risk score based on 5-mRNA signature significantly associates with prognosis. Distribution of each patient risk score (A, E, I, M) and individual mRNA expression profiles (D, H, L, P); the Kaplan–Meier survival analysis between high- and low-risk groups (B, F, J, N); time dependent ROC curves at 1, 2, and 3 years (C, G, K, O) in training set, internal test set, the entire TCGA set, and external test set. ROC = receiver operating characteristic, TCGA = The Cancer Genome Atlas.

Validation of the signature

To further assess the prognostic power of the identified mRNA signature, we determined its prognostic ability in internal TCGA and external datasets. Risk scores were evaluated using the training set formula. The distribution of the 5-mRNA risk score, patient survival, and mRNA expression were obtained in the internal test set (Fig. 2E and H), the entire TCGA set (Fig. 2I and L), and external test set (Fig. 2M and P). Kaplan–Meier curves demonstrated that patients in the low-risk group had improved outcomes compared to high-risk patients (Fig. 2F, J, and N). Comparable data were achieved through ROC analysis. The AUCs of the signature at 3-year survival time were 0.849, 0.924, and 0.747 in the internal test set (Fig. 2G), the entire TCGA set (Fig. 2K), and the external test set (Fig. 2O), respectively. These data were consistent with the training set data and validated the mRNA signature as a reliable predictor for overall survival in patients with EAC.

To investigate if the prognostic signature can be used independently of other risk factors, we performed univariable and multivariable Cox regression analyses across the TCGA dataset considering age, sex, height, weight, alcohol consumption, Barrett esophagitis, and TNM stage (Supplemental Table S2). As shown in the Figure 3, the mRNA signature was most significantly related to overall survival (hazard ratio, 4.32; 95% confidence interval, 1.91–9.79; P = .00046). In addition, the sensitivity and specificity of the mRNA signature performed numerically better than the TNM stage for 1-, 2-, and 3-year prognostic evaluation of EAC (Supplemental Fig. S3). When patients were stratified by clinicopathological factors, the mRNA signature represented a statistically significant prognostic model (Fig. 4). These findings demonstrate that high-risk scores possessed a strong association with a poor prognosis.

Figure 3

Forest plot summary of univariate and multivariate Cox regression analysis of overall survival in the whole TCGA set. CI = confidence interval, TCGA = The Cancer Genome Atlas.

Figure 4

Survival curves for EAC patients in different subgroup in the entire TCGA set. Kaplan–Meier analysis (low risk vs high risk) in subgroups stratified by age (A, B), sex (C, D), height (E, F), weight (G, H), alcohol consumption (I, J), Barrett esophagitis (K, L), TNM stage (M, N), and race (O). EAC = esophageal adenocarcinoma, TCGA = The Cancer Genome Atlas, TNM = tumor–node–metastasis.

Clinical and immune relevance of the mRNA signature

We next evaluated the correlation between the mRNA signature and clinical parameters (tumor size, node status, metastasis, and TNM stage) of EAC in TCGA set. Those with a positive lymph node status possessed higher risk scores (Fig. 5A). Stage III+IV patients also showed higher risk scores upon comparison to stage I+II patients (Fig. 5B).

Figure 5

Clinical- and tumor immunity relevance of the mRNA signature. The distribution of the risk score in different status of node and TNM stage in the whole TCGA set (A, B). The B cell infiltration level in high risk and low risk group in the entire TCGA set (C). The abundances of B cells were estimated using the TIMER algorithm (https://cistrome.shinyapps.io/timer/). The correlation of SINHCAF expression with B cell infiltration level in EAC patients (D). The Kaplan–Meier survival curves of different B cell infiltration level (E). ∗P < .05, ∗∗P < .01, ∗∗∗P < .001, ∗∗∗∗P < .0001. EAC = esophageal adenocarcinoma, TCGA = The Cancer Genome Atlas, TIMER = the tumor immune estimation resource, TNM = tumor–node–metastasis.

Using the TIMER analysis, the low-risk group was significantly related to higher levels of B cell infiltration (Fig. 5C). The expression of SLC25A45 negatively correlated with the infiltration of B cell (P < .0001; Fig. 5D). High levels of B cell infiltration prolonged survival in the EAC patients (P = .078; Fig. 5E).

Identification of biological processes associated with the mRNA signature

To infer the potential mechanisms of the mRNA signature, GSEA was applied to the TCGA and GSE72874 datasets. Gene sets based on KEGG and REACTOME pathways were associated with high-risk patients and identified in the TCGA datasets (Supplemental Table S3 and S4). The top 5 KEGG and REACTOME pathways in the TCGA are shown in Figure 6A and B. Meanwhile, GSEA was also performed in the GSE72874 dataset, and the results associated with KEGG and REACTOME pathways were presented in Supplemental Tables S5 and S6. Two commonly enriched KEGG pathways were screened including cell cycle and PI3K–Akt signaling pathway. The enriched REACTOME pathways were involved in nucleotide metabolism (DNA replication, DNA synthesis, and repair) and cell cycle progression (G1/S phase, S phase, and G2/M phase). We constructed enrichment maps and organized the enriched terms into networks with overlapping gene sets (Supplemental Fig. S4).

Figure 6

Functional analysis depicted the potential biological mechanism associated with the mRNA signature. The top 5 enriched KEGG pathways in patients with high risk (A). The top 5 enriched REACTOME pathways in patients with high risk (B). Significant GO terms enriched by the mRNA signature (C, D, E). GO = Gene Ontology, GSEA = Gene Set Enrichment Analysis, KEGG = Kyoto Encyclopedia of Genes and Genomes.

As shown in Figure 6C to E, the mRNA signature correlated with chromosomes in the CC, mediated MFs such as DNA-dependent ATPase activity and catalytic activity, and regulated DNA replication, chromosome segregation, and RNA splicing of the BP. These corresponding biological functions may contribute to poor prognosis of EAC patients with high-risk score.

Discussion

As a highly malignant neoplasm, the incidence of EAC is increasing and has surpassed ESCC in some areas of North America.[ TNM staging is a conventional and effective tool based on anatomical information. It helps us to improve current empirical treatment decisions and predict prognosis in patients with EAC. It is however unable to achieve adequate assessments of disease outcome in EAC. Molecular prognostic biomarkers such as mRNAs, microRNAs, and lncRNAs can supplement or substitute TNM staging in EAC.[ Gene expression and relevant clinical data of EAC available in TCGA and GEO facilitate the establishment of novel molecular signatures related to prognosis. Here, we employed TCGA and GEO datasets to build prognostic mRNA signatures for EAC.

We screened and selected significant DEMs associated with survival in the training set. Based on LASSO, we constructed a 5-mRNA signature in the training set and further validated its accuracy in internal and external test sets. As a supplement of TNM staging, the 5-mRNA signature could guide the stratification of patients with EAC and aid in decision making for individualized treatments.

SLC26A9, SINHCAF, MICB, KRT19, and MT1X formed the 5-mRNA signature that could predict the prognosis of EAC in the independent datasets. SINHCAF, also known as FAM60A, is a component of SIN3/HDAC deacetylase complex.[ In esophageal carcinoma, FAM60A was reported as a driver gene with a significant correlation with prognosis. FAM60A silencing could inhibit the proliferation, migration, and invasion of cells, increase apoptosis and arrest cells in the G2/M phase.[ FAM60A is overexpressed in gastric cancer tissues compared to adjacent tissue, and its overexpression enhances the gastric cancer cell proliferation through the PI3K/AKT pathway.[ In lung cancer and liver cancer cells, FAM60A acts as a tumor suppressor,[ which is not consistent with present study.

Major histocompatibility complex (MHC) class I chain-related protein B (MICB) is a ligand of NKG2D receptors on NK cells and T cells.[ When MICB is expressed on the cell surface, engagement between MICB and NKG2D actives the tumor killing effects of NK cells and T cells.[ In hepatocellular carcinoma and cervical cancer, patients with high expression of MICB showed an improved outcome.[ However, when MICB is retained in the cytoplasm, higher rates of immune evasion of cancer cells occurs, decreasing patient survival.[ The function of MICB markedly differs depending on the cancer type. Li et al[ reported that high levels of MICB are linked to poor prognosis in ovarian cancer. In ESCA, the expression of MICB was upregulated and associated with the histological grade of ESCA.[

The function of KRT19 in a range of cancer types differs to its role in EAC based on our findings. KRT19 promoted cancer cell proliferation and invasion and was identified as a marker of poor prognosis for hepatocellular carcinoma, pancreatic ductal adenocarcinoma, and breast cancer.[ Liu et al[ found that MT1X could serve as a favorable prognostic marker and inhibit the progression and metastasis of hepatocellular carcinoma. In oral squamous cell carcinoma, the upregulation of MT1X was related with a lack of metastasis which indicated that MT1X may aid prediction of prognosis.[ However, few cancer-related studies have been performed regarding the function of SLC26A9, which warrants further investigation.

In this study, we found that B cell infiltration was significantly higher in low-risk groups and suggestive of an improved outcome. T cells are therapeutic targets with the immune checkpoint inhibitors in ESCA, but the efficacy of this therapy is variable.[ B cells are major effector cells of humoral immunity and mediate immune responses. They prevent cancer progression through immunoglobulins secretion, T cell responses enhancement and cancer cell killing.[ To balance immunotherapy and adverse events, further stratification based on B cell infiltration or risk scores are necessary to make informed treatment decisions. Further experimental studies are similarly required to decipher the function of B cells during cancer progression.

To our knowledge, this mRNA signature associated with prognosis built using the TCGA and GEO datasets have not been reported previously. Inevitably, there are some limitations in the current study. First, although 2 datasets were used to construct and validate the signature, the sample size was small. This signature requires further validated in large prospective clinical trials. In addition, Dong et al[ demonstrated the function of SINHCAF in ESCA, but no experimental studies on the other genes of the signature in ESCA have been reported. To better understand the potential mechanisms behind the signature, in vitro and in vivo experiments should be performed.

In conclusion, despite the limitations described, the 5-mRNA signature that we constructed through a comprehensive analysis of the mRNA profiles of EAC may serve as a novel and reliable tool to aid prognosis predictions and tailored treatments for EAC. This study also provides significant information for molecular research and clinical treatment in EAC.

Acknowledgments

The authors sincerely thank the TCGA and GEO datasets for offering access to their resources.

Author contributions

Conceptualization: Tian Lan, Hua Luo.

Supervision: Hua Luo.

Validation: Yunyan Lu.

Visualization: Yunyan Lu.

Writing – original draft: Tian Lan.

Writing – review & editing: Tian Lan, Weiguo Liu, Hua Luo.