Identification the prognostic value of immune gene signature and infiltrating immune cells of esophageal cancer patients

doi:10.21203/rs.3.rs-20941/v1

Download PDF

Research

Identification the prognostic value of immune gene signature and infiltrating immune cells of esophageal cancer patients

https://doi.org/10.21203/rs.3.rs-20941/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 01 Oct, 2020

Read the published version in International Immunopharmacology →

Version 1

posted

You are reading this latest preprint version

Background Esophageal cancer (ESCA) is one of the deadliest solid malignancies with worse survival in the world. The poor prognosis of ESCA is not only related to malignant cells, but also affected by the microenvironment. We aimed to establish prognostic signature consisting of immune genes to predict the survival outcome of patients and estimate the prognosis value of infiltrating immune cells in tumor microenvironment (TME).

Methods Based on integrated analysis of gene expression profiling and immune gene database, differentially immune-related genes were filtered out. Then, stepwise Cox regression analysis was applied to identify survival related immune genes and construct prognosis signature. Functional enrichment analysis was performed to explore biology function. Kaplan-Meier (K-M) curves and receiver operating characteristic (ROC) curves were performed to validate the predictive effect of predictive signature. We also verified the clinical value of prognostic signature under the influence of different clinical parameters. For deeper analysis, we evaluated the correlation between prognosis signature and infiltrating immune cells by Tumor Immune Estimation Resource (TIMER) and CIBERSORT.

Results Finally, we identified 303 differentially immune genes as candidate and constructed immune prognosis signature composed of six immune genes. Furthermore, we observed that the prognosis signature was enriched in cytokine-mediated signaling pathway, lymphocyte activation, immune effector process, cancer pathway, NF-kappa B signaling pathway. K-M survival curves showed that the prognosis signature indeed have good predictive ability in entire ESCA set ( P =0.003), validation set 1 ( P =0.008) and validation set 2 ( P =0.036). The area under the curve (AUC) of ROC curves validated the predictive accuracy of immune signature in three cohorts (AUC=0.757, 0.800 and 0.701), respectively. In addition, we identified the prognosis value of infiltrating-immune cells including activated memory CD4 T cells, T cells follicular helper cells and monocytes and provided a landscape of TME.

Conclusions The results indicated that immune prognosis signature can be a novel biomarker to predict survival outcome, which can provide new targets for immunotherapy and individualized therapies in ESCA and open up a new prospect for improving the prognosis of ESCA patients in the era of immunotherapy.

Translational Medicine

immune related genes

infiltrating immune cells

prognosis signature

esophageal cancer

tumor microenvironment

Esophageal cancer (ESCA) is one of the most deadly cancers worldwide which incidence rate is increasing by annually [1]. Although the treatment of ESCA has been improved, the 5-year survival rate is about 15% which lead to worse prognosis and higher mortality [2]. For ESCA patients, due to tumor heterogeneity, even patients with the same clinical conditions will have different results after the same treatment [3]. The prognosis of patients can’t achieve a good judgment and missing the best time and method for treatment. Accurate evaluation the prognosis of patients can provide individualized treatment for each patient and improve their survival rate [4]. Therefore, it is imperative and significant to identify prognosis signature for survival prediction of ESCA patients. According to prognosis biomarkers, we can choose the personal therapy and predict the survival outcome of patients.

Biomarkers based on genes have gradually become widespread application and more cost-effective way to predict prognosis of ESCA patients [5]. Gene expression profiling have revealed the presence of inflammation which have huge value for the study of prognosis markers related to tumor survival [6]. Notably, the deterioration of prognosis in ESCA will increase along with genes changes. Similarly, the whole genome sequencing research also pointed out that the gene profiling of ESCA patients is different from that of normal people [7]. Researchers have realized that gene expression profile plays an extremely irreplaceable role in the analysis of cancer development, which can make poor survival rate in ESCA [8]. However, the value of single gene biomarker in predicting the prognosis is limited because of the tumor heterogeneity, especially caused by the alteration of genes. The application of integrated biomarkers in esophageal squamous cell carcinoma (ESCC), one of the subtypes of ESCA, has been reported in succession. For instance, DNA methylation related five-gene signature was identified in ESCC [9]. Mao and colleagues established seven-lncRNA signature for survival prediction in ESCC [10]. Therefore, gene signatures consist of multiple genes are becoming a better choice and it’s necessary to identify effective robust combined biomarkers that can indicate prognosis in ESCA patient [11].

At present, immunotherapy and targeted therapy are wide application and become an important way to improve the survival rate of ESCA patients. Different from the traditional treatment, the target of immunotherapy is not cancer cells, but some immune checkpoint including cytotoxic T-lymphocyte associated antigen 4 (CTLA-4) and programmed death 1 (PD-1) [12, 13]. However, some patients are not sensitive to inhibitors of immune checkpoint. The great challenge of tumor immunotherapy is to identify prognosis models to determine the therapeutic effect and guide the treatment of diseases [14]. Therefore, it is of great significance for the individualized treatment of ESCA to find effective biomarkers that can predict the sensitivity of patients to immunotherapy. Studies reported that immune related genes are not only correlation with response of immunotherapy, but also related to the prognosis of patients. For instance, the prognostic signature of seven immune genes was identified in clear cell renal cell carcinoma [15]. An immune gene-set can predict the prognosis of patients with ovarian cancer [16]. Guo and colleagues established immune signature to predict survival of lung adenocarcinoma [17]. Although many studies have explored the relationship between immune related genes and prognosis of patients, few studies to establish prognosis signature of immune genes in ESCA based on expression profile data. Therefore, a reliable prognosis model composed of immune genes must be established for ESCA and reveal the value of clinical application.

Besides of genomic features, tumor microenvironment (TME) represents promising candidates for predictive and prognostic biomarkers [18]. TME plays a critical role in cancer initiation, progression, therapeutic response and clinical reaction [19]. In the wake of the deeper research, the type and density of infiltrated immune cells in TME reported can regulate tumor cells to induce host immune response and significantly influence the development of cancer [20]. Tumor-infiltrating immune cells including B cells, macrophages (M0, M1 or M2), effector T cells, natural killer cells, mast cells, naive and memory lymphocytes and dendritic cells may be found around tumor cells in TME. Different infiltrating immune cells in ESCA patients will interact with tumor cells in various ways to resist treatment. Therefore, the integrated analysis of prognosis signature and infiltrating-immune cells is of great significance to predict the response of immunotherapy and the prognosis of patients.

In this study, we integrated gene expression data set and immune gene database to verify survival related immune genes and construct immune prognostic signature in ESCA. Meanwhile, we validated the clinical value of prognostic signature under the influence of different clinical parameters. Given that tumor immune microenvironment related to prognostic features may be involved in the development of new biomarkers, we discussed the correlation between identified prognosis signature and immune cell infiltration. Then, we validated the prognostic effect of immune cell subtypes between risk groups defined by immune risk model. The results indicated that immune prognosis model related to TME can be a novel biomarker to predict survive, which will open up a new prospect for improving the prognosis of ESCA patients in the era of immunotherapy.

Data acquisition

The transcriptome profiling of RNA expression [FPKM] and clinical data of ESCA patients consisting of 160 tumor samples and 11 normal samples were collected from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/). The TCGA data provide a valuable source to analyze complex cancer genomics and clinical parameters. We removed one patient whose transcriptomic data and clinical data are not completed. Thus, the entire TCGA data (n = 159) was as the training set for further analyses. The clinical information of patients was summarized in Table 1. Furthermore, 2498 immune-related genes were obtained for further immunology research from Immunology Database and Analysis Portal (ImmPort, https://www.immport.org/home) [21]. It includes 17 immune categories such as T-cell receptor signaling pathway, interleukins, B-cell receptor signaling pathway, cytokine, tumor necrosis factor (TNF) family receptors and so on. Therefore, TCGA database and ImmPort database were used to develop prognosis gene signature.

Table 1

Summary of clinical characteristics of ESCA patients in three cohorts
Characteristic	Patients in entire TCGA set (n = 159), n (%)	Patients in validation set 1 (n = 80), n (%)
Age(years)
≤ 60	81 (50.94%)	33 (41.25%)	48 (60.76%)
> 60	78 (49.06%)	47 (58.75%)	31 (39.24%)
Gender
Female	23 (14.47%)	16 (20.00%)	7 (8.86%)
Male	136 (85.53%)	64 (80.00%)	72 (91.14%)
Histological type
EAD	79 (49.69%)	40 (57.14%)	39 (49.37%)
ESCC	80 (50.31%)	30 (42.86%)	40 (50.63%)
Barretts esophagus
NO	104 (80.62%)	57 (86.36%)	47 (74.60%)
YES	25 (19.38%)	9 (13.64%)	16 (25.40%)
Vital status
Alive	96 (60.38%)	53 (66.25%)	43 (54.43%)
Dead	63 (39.62%)	27 (33.75%)	36 (45.57%)
Neoplasm Cancer Status
Tumor Free	91 (61.10%)	46 (61.33%)	45 (60.81%)
With Tumor	58 (38.90%)	29 (38.67%)	29 (39.19%)
Pathologic stage
Stage I-II	87 (56.13%)	43 (55.84%)	44 (56.41%)
Stage III-IV	68 (43.87%)	34 (44.16%)	34 (43.59%)
Alcohol history
NO	46 (29.49%)	23 (29.11%)	23 (29.87%)
YES	110 (70.51%)	56 (70.89%)	54 (70.13%)
T Classification
T1-T2	66 (45.14%)	36 (46.15%)	30 (37.97%)
T3-T4	81 (54.86%)	42 (53.85%)	49 (62.03%)
N Classification
N0	64 (43.54%)	30 (41.10%)	34 (45.95%)
N1-N3	83 (56.46%)	43 (58.90%)	40 (54.05%)
M Classification
M0	126 (89.36%)	66 (91.67%)	60 (86.96%)
M1	15 (10.64%)	6 (8.33%)	9 (13.04%)
New event
NO	89 (55.97%)	49 (61.25%)	40 (50.63%)
YES	70 (44.03%)	31 (38.75%)	39 (49.37%)
Tumor Central Location
Proximal and Mid	47 (29.75%)	26 (32.50%)	21 (26.92%)
Distal	111 (70.25%)	54 (67.50%)	57 (73.08%)
Neoplasm Histologic Grade
G1	16 (12.90%)	10 (17.24%)	6 (9.09%)
G2	65 (52.42%)	30 (51.73%)	35 (53.03%)
G3	43 (34.68%)	18 (31.03%)	25 (37.88%)
Lymph node metastasis
NO	83 (65.87%)	43 (68.25%)	40 (63.49%)
YES	43 (34.13%)	20 (31.75%)	23 (36.51%)
Radiation therapy
YES	16 (14.68%)	11 (18.18%)	5 (7.81%)
NO	93 (85.32%)	54 (81.82%)	59 (92.19%)

Identification of differentially immune genes

At first, we obtained differentially expression genes from TCGA data. In order to identify differentially genes between tumor tissues and normal tissues, the expression differences were characterized by absolute log 2 fold change (FC) > 1 and adjusted P-value < 0.05. Next, intersection differentially expression genes with immune-related genes from ImmPort database to screen differentially expressed immune-related genes. Finally, differentially immune-related genes were filtered out totally. These selected genes as initial candidate were used to establish the immune-related risk signature for next step. In addition, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of differentially immune genes were performed by Metascape (http://metascape.org). This website can obtain biological pathways through independent and orthogonal experiments on data sets of more than 40 knowledgebase [22]. In general, P ≤ 0.05 were represented significantly enriched pathways.

Construction of the immune-related risk score system

We further investigated these selected genes for constructing immune related prognosis signature. Univariate Cox regression analysis was applied to identify possible prognosis immune genes. Multivariate linear Cox hazards regression was performed to construct the best prognostic model in ESCA patients by survival package in R software. Finally, the prognostic signature (six-gene), risk score, was calculated based on gene expression and Cox regression coefficient as followed:

In the formula, n, Expi, and βi represented the number of hub genes, gene expression level and regression coefficient value, respectively. Accordingly, patients were ranked with risk score into high and low risk group by using the median value as the cutoff point. The distribution of the survival status was validated based on risk score levels. The higher risk score, the poor prognosis ESCA patients have.

Assessment the prognosis value of immune signature

The 159 patients were randomly divided into two sets: validation set 1 (n = 80) and validation set 2 (n = 79) (Table 1). The prognostic signature was identified in entire TCGA data set and the prognosis performance of the model was validated in both three sets. Given the effect of immune signature on survival outcomes in ESCA, we then verified the prognosis value of the six-gene model. K-M survival curves with Log-Rank test were applied to compare the OS effect of prognosis signature in high and low risk groups. The predictive accuracy of prognosis signature was determined by the time-dependent receiver operating characteristic (ROC) curve [23]. The area under the curve (AUC) of ROC curve was calculated to evaluate the accuracy of the immune signature by survival ROC package in R and an AUC > 0.70 was considered to have good predictive value. Likewise, we used uni- and multivariate linear Cox regression analysis to investigate the influence of clinical parameters in the prognosis signature by survival package in R and IBM SPSS 25.0 program.

Exploration the prognosis value of infiltrating immune cells in TME

Given the key effect of infiltrating immune cells in TME, we analyzed the relationship between immune genes and infiltrating immune cells by Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/). The website of TIMER can be used to analyze the fraction of six immune infiltrating cells including B cells, CD8⁺ T cells, CD4⁺ T cells, neutrophils, macrophages and dendritic cells in TME, and the correlation between genes and different immune cells. Furthermore, CIBERSORT algorithm was performed to identify the fraction of infiltrating immune cells (n = 22) in patients with ESCA. CIBERSORT (http://cibersort.stanford.edu/) is a computational framework which can estimate the proportions of 22 distinct immune cell types with gene expression profiles from TCGA-ESCA [24]. In final, a total number of 98 patients with CIBERSORT P value > 0.05 and OS were included. The intersection between 159 patients and 98 patients with risk score showed 98 patients were analyzed for further study. Based on risk score of prognosis signature, 98 patients were divided into high and low risk groups by using the median cutoff value. In order to validate prognosis value of infiltrating immune cells, K-M curves and log-rank test were performed in ESCA patients. Finally, we identified immune cells which have the function to predict patients’ survival status. We also investigated the potential small molecule drugs for immunotherapy of ESCA patients by drug-bank (https://www.drugbank.ca/drugs). Based on the public websites or databases, we could discover immune cells related to survival and reveal the relationship between immune prognosis signature and immune infiltration cells in TME. Statistically significant P-value was set as ≤ 0.05.

Statistical analysis

In the present, Student t-test (for equal variances) was performed to statistical comparison. GraphPad Prism 8.0 software was used to plot K-M survival curves. The analysis of Cox regression was conducted by IBM SPSS 25.0 program. The heatmaps and pheatmaps were generated by applying R package. All statistical analysis was conducted to determine independent prognostic factors which can predict patients’ survival status by the R package (R software version 3.5.2). P-value (two-sided test) less than 0.05 was considered significant with the purpose of ensuring the reliability of the results.

Identification of differentially immune genes in patients of ESCA

In this study, we obtained 4102 differentially expression genes from TCGA-ESCA data in totally (Additional file 1: Figure S1a-b). Then, we intersected these genes with 2498 immune-related genes from ImmPort database, which was performed by Venn diagram (Additional file 1: Figure S1c). In the end, 303 differentially immune-related genes were filtered out for further analysis (|log 2 FC| > 1, P-value < 0.05, Additional file 1: Figure S1d-e). In order to exploration of the potential biological significant of the differentially expressed immune genes, GO categories (top 20, Additional file 2: Figure S2a, Additional file 3: Table S2) and KEGG pathway (top 20, Additional file 2: Figure S2b, Additional file 4: Table S3) were applied by Metascape. In GO analysis, we detected that these differentially immune genes are mainly enriched in immune related pathways including receptor ligand activity, cytokine-mediated signaling pathway, response to tumor necrosis factor, lymphocyte activation, response to interferon-gamma, immune effector process and cell proliferation. The significant KEGG pathways were also enriched in immune related signal pathway like chemokine and cytokine signaling pathway, natural killer cell mediated cytotoxicity, cancer pathway, NF-kappa B signaling pathway, T cell or Th17 cell signaling pathway and PI3K-Akt signaling pathway. Meanwhile, we represented these pathways in a network way, which is more conducive to observe and understand the interaction of these biological processes (Additional file 5: Figure S3). Taken together, these signal pathways represented that these immune genes (n = 303) play a crucial role in immune and cancer responses.

Construction of the immune-related prognosis signature in the entire TCGA set

For construction of prognosis signature, the entire ESCA patients (n = 159) as training cohort. Firstly, survival related immune genes were selected by performing univariate Cox hazards regression analysis. We thus obtained survival genes as candidate for further construction of prognostic signature in ESCA patients. Furthermore, multivariate linear Cox hazards regression analysis were applied to select independent prognostic factors of expressed immune genes. Finally, we established prognosis signature consisting of six immune-related genes (HSPA6, S100A12, FABP3, DKK1, OSM, and NR2F2, Table 2). Meanwhile, risk score was calculated based on gene expression and Cox regression coefficient as followed:

Risk score = 0.008 * Exp of HSPA6 + 0.004 * Exp of S100A12+0.042 * Exp of FABP3 + 0.010 * Exp of DKK1 + 0.271 * Exp of OSM + 0.012 * Exp of NR2F2

Table 2

The detailed information of identified six genes significantly related to overall survival in ESCA patients
ID	Coef(β_i)	HR	95%CI of HR	P-value
HSPA6	0.008	1.008	1.004–1.013	0.000
S100A12	0.004	1.004	1.002–1.005	0.000
FABP3	0.042	1.043	0.994–1.093	0.085
DKK1	0.010	1.010	1.000-1.020	0.049
OSM	0.271	1.312	1.119–1.537	0.001
NR2F2	0.012	1.012	1.001–1.023	0.032

Then, entire ESCA patients were randomly divided into two sets: validation set 1 (n = 80) and validation set 2 (n = 79) to validate the risk model. Based on median of risk value, patients were separated into low-risk and high-risk groups in entire TCGA set (n = 159, Fig. 2a), TCGA test set (n = 80, Fig. 2b) and TCGA validation set (n = 79, Fig. 2c), respectively. We found that patients with increasing risk score have shorter survival time both in three cohorts. These results showed that six immune-gene signature (risk score system) may have more prognosis function. Given the significance of prognosis model to clinical application, we performed univariate and multivariate survival tests to prove the relationship between risk model and different clinical pathological parameters including pathological stage-N, pathological stage-M, cancer stage, cancer status, lymph node metastasis, radiation therapy and new event (Table 3). Univariate survival analyses showed that risk score and some clinical parameters (pathological stage-N, pathological stage-M, cancer stage, cancer status and lymph node metastasis) can be prognosis biomarker in ESCA patients (n = 159, Fig. 2d). Further multivariate survival analysis confirmed that risk score (P = 0.003) and cancer status (P < 0.001) were independent prognostic indicators for ESCA patients (Fig. 2e). Although other clinical indicators in multivariate survival analysis were less powerful, they still have potential value in clinical application. It is notable that risk score risk model also has prognostic significance under the influence of clinical factors.

Table 3

Univariable and multivariable Cox regression analyses for risk score and different clinical pathological parameters in ESCA
		Univariable analyses			Multivariable analyses
Clinical feature	Number	HR	95%CI of HR	P-value	HR	95%CI of HR	P-value
Risk score(low risk/ high risk)	81/78	2.198	1.307–3.697	0.003	3.203	1.472–6.967	0.003
Pathological stage-N(N0/N1-N3)	64/83	3.018	1.611–5.651	0.001	2.774	0.958–8.032	0.060
Pathological stage-M(M0/M1)	126/15	3.159	1.656–6.027	< 0.001	0.468	0.158–1.391	0.172
Cancer stage(stage I-II/stage III-IV)	87/68	2.759	1.626–4.679	< 0.001	1.322	0.498–3.509	0.575
Cancer status(tumor free/with tumor)	91/58	3.170	1.817–5.531	< 0.001	10.055	2.939–34.403	< 0.001
Lymph node metastasis(no/yes)	83/43	1.846	1.061–3.211	0.030	0.513	0.184–1.427	0.201
Radiation therapy(yes/no)	16/113	0.797	0.374–1.698	0.557	2.884	0.835–9.966	0.094
New event(no/yes)	89/70	1.261	0.764–2.082	0.364	0.222	0.085–0.584	0.002

Validation of the predictive performance of immune prognosis signature

Based on uni- and multivariate Cox hazards regression analysis, we observed that prognosis signature can be independent prognostic biomarker to indicate patients’ survival outcome. In order to validate the predictive capability of prognostic signature, survival curves were performed by K-M survival analysis in three cohorts. The time-dependent ROC curves were applied to assess the accuracy of predictive function of this prognosis model in three cohorts. In entire TCGA data set, the K-M survival curve showed that the prognostic signature indeed can well distinguish patients into high or low survival rate (P = 0.0034, Fig. 3a). ROC curve (AUC) for predicting patients survival confirmed that the identified prognostic signature has the robust efficiency to predict the OS for ESCA patients (AUC = 0.757, Fig. 3b). In validation set 1, the K-M survival curve (P = 0.0081, Fig. 3c) and ROC curve (AUC = 0.800, Fig. 3d) also showed that the prognosis signature indeed have good predictive ability in ESCA patients. In validation set 2, the K-M survival curve also showed evident gaps between low-risk and high-risk patients (P = 0.0363, Fig. 3e) and the ROC curve (AUC = 0.701, Fig. 3f) again validated. We observed that the median survival time (MST) of patients in high risk group was lower than other group in both three cohorts. These results showed that the prediction model can be good prognosis indicator in patients with ESCA.

Independent prognostic value of six immune-gene signature and clinical parameters

Based on the previous stepwise Cox regression analysis, we know that clinical parameters were effective prognosis predictors in patients with ESCA. Thus, we performed the K-M survival curves to validate the prognosis difference of clinical parameters including pathological stage-N (P = 0.0001, Fig. 4a), pathological stage-M (P = 0.0002, Fig. 4b), cancer stage (P < 0.0001, Fig. 4c), cancer status (P < 0.0001, Fig. 4d) and lymph node metastasis (P = 0.0196, Fig. 4e) in two (high/low) risk groups. Obviously, clinical parameters can distinguish survival difference between different subtypes. The difference of MST between two subtypes of the same clinical factor was significant. Patients who were in pathological stage-N0, in pathological stage-M0, in stage I-II, tumor free (without tumor) and without lymph node metastasis had better prognosis compared with other subtypes. These results further validated that these clinical indicators indeed have good prognosis value in patients of ESCA.

Considering the clinical value of prognostic model, stratified survival analysis was performed to validate whether the prognosis model can distinguish survival difference between different subtypes. For clinical parameters of pathological stage-N (N0 or N1-N3, Fig. 5a-b), pathological stage-M (M0 or M1, Fig. 5c-d), cancer stage (stage I-II or stage II-IV, Fig. 5e-f), cancer status (tumor free or with tumor, Fig. 5g-h) and lymph node metastasis (No or Yes, Fig. 5i-j), we found that the MST of patients with high risk was lower than other group. Risk signature based on the integrated six immune genes can be confirmed as independent prognostic indicator and indeed significant in clinical application in ESCA.

Revealing the relationship between prognostic signature and tumor-infiltrating immune cells in TME

Study reported, infiltration of immune cells in tumor microenvironment was accompanied by cancer initiation and progression [25]. Further investigation indicated that the presence of infiltrating immune cells can be used as biomarker for immunotherapy response [26]. Therefore, we evaluated the correlation between prognosis signature we identified and six types of immune cells. Based on the database of TIMER, we discovered that six immune genes from prognosis signature have good relationship with these infiltrating immune cells (Additional file 6: Figure S4). In order to estimate more immune cells, we applied CIBERSORT algorithm [27] to verify the proportions of 22 immune cells ESCA. Firstly, we assessed the abundance of distinct immune cells in each ESCA patient (Fig. 6a) and two (low/high) risk groups (Fig. 6b) in bar charts. Then, we investigated the correlation of each cell type which showed that B cells naive were highly correlated with plasma cells in ESCA patients (Fig. 6c). In order to observe the prognosis value of immune cells, K-M survival curves were performed in patients. Survival curves confirmed that activated memory CD4 T cells (P = 0.0195, Fig. 6d), T cells follicular helper (P = 0.0380, Fig. 6e) and monocytes (P = 0.0418, Fig. 6f) have prognostic benefit for ESCA patients. All above results show that the infiltrating immune cells in TME were associated with prognosis. In addition, based on drug targeted analysis of prognosis signature, we found that Olopatadine and Amlexanox can be potential drugs of immunotherapy in ESCA patients by drug-bank database. All of them can target gene of S100A12 (one of six hub genes in prognosis signature). Studies reported the inhibition of inflammation for these drugs [28, 29]. We deeply hope these small molecule drugs have clinical value and applied for patients to improve prognosis. Therefore, what we identified in the study could provide prognosis biomarker based on immune genes for ESCA patients and eventually apply to personalized immune-targeted therapy.

ESCA is an aggressive disease with highly malignant and poor prognosis [30]. In recent years, novel biomarkers established to predict the prognosis of cancer patients have been emerging, especially immune related gene signatures. Although immunoassays have been performed in a variety of cancers, related analysis in ESCA is rare. We should carry on the prognosis analysis of immune genes in ESCA patients to explore clinical significance and potential molecular mechanism. In our study, we developed a prognostic signature consisting of multiple immune genes to predict the survival outcome by integrated analysis of various databases. Finally, we construct a prediction model with six genes (HSPA6, S100A12, FABP3, DKK1, OSM, and NR2F2), which can predict the prognosis powerfully. In addition, we detected immune cells infiltration and evaluate potential clinical results.

Notably, all six hub genes in signature have prognosis value in the development of cancers. Study reported that HSPA6 (heat shock proteins, member of HSP70 family) was related to suppress apoptosis, cancer and participate in TME, which was identified as potential targets in ovarian cancer [31]. Overexpression of HSPA6 played a critical role in the recurrence of hepatocellular carcinoma (HCC), which was validated as a potential biomarker in HCC [32, 33]. S100A12 belongs to calcium-binding S100 group proteins, which was necessary in cancer development. S100A12 was related to the proliferation and invasion of human papillary thyroid cancer (PTC), and may be an effective target for PTC treatment [34]. Meanwhile, S100A12 can be important diagnostic biomarker in breast cancer [35]. It has been reported that S100A12 have prognosis value in HCC [36]. FABP3 is fatty acid binding proteins, which involved in tumor genesis [37]. Tang and colleagues have reported that high expression of FABP3 can indicate worse survival in patients of non-small cell lung cancer [38]. Studies also reported that FABP3 can be prognosis biomarker in breast cancer [39] and ovarian cancer [40]. DKK1 (Dickkopf-related protein 1) was associated with prognosis of breast cancer [41]. It has been demonstrated that DKK1 can be prognosis factor in pancreatic ductal adenocarcinoma [42]. OSM (oncostatin M) was involved in tumor development and insulin resistance [43]. It has been reported that OSM as a risk factor promote breast cancer metastasis [44]. Studies indicated that the expression of NR2F2 (Nuclear receptor subfamily 2) are associated with survival outcome in gastric cancer [45] and breast cancer [46]. In general, we found that these genes are related to cancer or prognosis by reviewing the existing reports. The prognosis signature based on six genes is very reliable.

In this study, we established a TME-related prognosis biomarker to predict survival ending for ESCA patients based on immune related genes by integrated analysis of databases. Afterwards, we performed univariate and multivariate Cox regression analysis to identify survival related immune genes. Meanwhile, patients were separated to two (low/high) risk groups by risk score system. Finally, we developed six-gene biomarker as a novel prognosis model and analyzed their ability to predict prognosis in different cohorts. The results have confirmed that the predictive model of immune genes could independently distinguish that high risky patients have worse survival compared to low risky group in ESCA. Given the importance of clinical factors, K-M survival curves of OS were applied in different clinical pathological parameters. As shown in results, clinical factors are associated with OS of patients. Importantly, the prognosis signature is also independent indicator and has good performance to predict the survival when adjustment clinical parameters including pathological stage-N (N0 or N1-N3), pathological stage-M (M0 or M1), cancer stage (stage I-II or stage II-IV), cancer status (tumor free or with tumor) and lymph node metastasis (No or Yes). In conclusion, the six-gene model can predict patients’ prognostic features and apply on ESCA patients for appropriate therapy to prolong their survival time.

The accelerated deterioration of cancer is not only related to malignant cells, but also affected by the microenvironment [47]. Tumor-infiltrating immune cells in TME play an important role in improving antitumor and immunotherapeutic effects. Therefore, we discussed the relationship between immune gene signature and infiltrating immune cells to reflect the TME of ESCA. Notably, six hub genes of signature have relationship with many immune cells including B cells, macrophage, dendritic cells and T cells by performed TIMER (only including six types of immune cells). Furthermore, we applied CIBERSORT to detected survival-related immune cells. In results, activated memory CD4 T cells, T cells follicular helper cells and monocytes can predict survival of patients. The importance of T cells and monocytes in cancer immunotherapy was indeed reported [48, 49]. Study reported that T cells follicular helper cells can indicate prognosis in BRCA patients [50]. Therefore, not only genomic signature but also infiltrating immune cells in TME can represent promising candidates for predictive and prognostic biomarkers.

In conclusion, we applied an integrated study and established a novel six immune gene signature for the survival prediction of ESCA patients. We verified that the prognosis signature has clinical implications in ESCA patients. In addition, the prognosis signature provided a landscape of TME and revealed OS related immune cells. Taken together, our study could contribute to elucidate the prognosis value of gene signature and provided new targets for immunotherapy and individualized therapies in ESCA.

ESCA: esophageal cancer

ESCC: esophageal squamous cell carcinoma

CTLA-4: cytotoxic T-lymphocyte associated antigen 4

FC: fold change

PD-1: programmed death 1

TNF: tumor necrosis factor

ImmPort: Immunology Database and Analysis Portal

TME: tumor microenvironment

TCGA: The Cancer Genome Atlas

OS: overall survival

ROC: receiver operating characteristic

K-M: Kaplan-Meier

MST: median survival time

GO: Gene Ontology

KEGG: Kyoto Encyclopedia of Genes and Genomes

TIMER: Tumor Immune Estimation Resource

Ethics approval and consent to participate

All data in this study were performed after local institution review boards approved and consent to participate.

Consent for publication

The participant has consented to the submission to the journal.

Conflicts of Interest

Authors declare no conflicts of interest.

Data Availability Statement

All data included in this study are available including The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/) and Immunology Database and Analysis Portal (ImmPort, https://www.immport.org/home).

Authors’ contributions

LW designed this study. LW, QW and MH collected and processed all data. LW, MZ and LZC analyzed all data. LW, ZNL and CYZ prepared tables and figures for this study. LW drafted and achieved the manuscript. MJW and LZ revised the manuscript. All authors read and approved the final manuscript.

Acknowledgements

The authors would like to acknowledge the Key Laboratory of Precision Diagnosis and Treatment of Gastrointestinal Tumors, Ministry of Education (China Medical University, Shenyang, China) for providing the space and equipment for conducting the experiments.

Funding

This work was supported by grants from National Natural Science Foundation of China (No. U1608281), Liaoning Revitalization Talents Program (No. XLYC1807201) and Shenyang S&T Projects (No. 19-109-4-09).

Author information

Affiliations

Department of Pharmacology, School of Pharmacy, China Medical University, Shenyang, Liaoning Province, China

Liaoning Key Laboratory of molecular targeted anti-tumor drug development and evaluation; Liaoning Cancer immune peptide drug Engineering Technology Research Center; Key Laboratory of Precision Diagnosis and Treatment of Gastrointestinal Tumors, Ministry of Education; China Medical University, Shenyang, Liaoning Province, China.

Corresponding authors

Correspondence to Minjie Wei and Lin Zhao

E-mail address: [email protected] (Minjie Wei)

E-mail address: [email protected]. (Lin Zhao)

Siegel RL, Miller KD, Jemal A: Cancer statistics, 2020. CA Cancer J Clin 2020, 70:7-30.
Domper AM, Ferrandez AA, Lanas AA: Esophageal cancer: Risk factors, screening and endoscopic treatment in Western and Eastern countries. World J Gastroenterol 2015, 21:7933-7943.
Lagergren J, Smyth E, Cunningham D, Lagergren P: Oesophageal cancer. Lancet 2017, 390:2383-2396.
J L, X L, C Z, C Z, H W: A signature of tumor immune microenvironment genes associated with the prognosis of non‑small cell lung cancer. Oncology reports 2020, 43:795-806.
He W, Yan Q, Fu L, Han Y: A five-gene signature to predict the overall survival time of patients with esophageal squamous cell carcinoma. Oncol Lett 2019, 18:1381-1387.
YQ L, F W, JJ L, YF L, X L, Z W, RC C: Gene Expression Profiling Stratifies IDH-Wildtype Glioblastoma With Distinct Prognoses. Frontiers in oncology 2019, 9:1433.
He W, Chen L, Yuan K, Zhou Q, Peng L, Han Y: Gene set enrichment analysis and meta-analysis to identify six key genes regulating and controlling the prognosis of esophageal squamous cell carcinoma. J Thorac Dis 2018, 10:5714-5726.
Zhang H, Zhong J, Tu Y, Liu B, Chen Z, Luo Y, Tang Y, Xiao F, Zhong J: Integrated Bioinformatics Analysis Identifies Hub Genes Associated with the Pathogenesis and Prognosis of Esophageal Squamous Cell Carcinoma. Biomed Res Int 2019, 2019:2615921.
Lu T, Chen D, Wang Y, Sun X, Li S, Miao S, Wo Y, Dong Y, Leng X, Du W, Jiao W: Identification of DNA methylation-driven genes in esophageal squamous cell carcinoma: a study based on The Cancer Genome Atlas. Cancer Cell Int 2019, 19:52.
Mao Y, Fu Z, Zhang Y, Dong L, Zhang Y, Zhang Q, Li X, Liu J: A seven-lncRNA signature predicts overall survival in esophageal squamous cell carcinoma. Sci Rep 2018, 8:8823.
Visser E, Franken IA, Brosens LA, Ruurda JP, van Hillegersberg R: Prognostic gene expression profiling in esophageal cancer: a systematic review. Oncotarget 2017, 8:5566-5577.
Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, Powderly JD, Carvajal RD, Sosman JA, Atkins MB, et al: Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med 2012, 366:2443-2454.
Pardoll DM: The blockade of immune checkpoints in cancer immunotherapy. Nat Rev Cancer 2012, 12:252-264.
Paucek RD, Baltimore D, Li G: The Cellular Immunotherapy Revolution: Arming the Immune System for Precision Therapy. Trends Immunol 2019, 40:292-309.
C S, J L, J W, X Z, D D, X Y, Y W: Development and validation of a prognostic immune-associated gene signature in clear cell renal cell carcinoma. International immunopharmacology 2020, 81:106274.
Shen S, Wang G, Zhang R, Zhao Y, Yu H, Wei Y, Chen F: Development and validation of an immune gene-set based Prognostic signature in ovarian cancer. EBioMedicine 2019, 40:318-326.
D G, M W, Z S, J Z: A new immune signature for survival prediction and immune checkpoint molecules in lung adenocarcinoma. Journal of translational medicine 2020, 18:123.
Chen Q, Liu G, Liu S, Su H, Wang Y, Li J, Luo C: Remodeling the Tumor Microenvironment with Emerging Nanotherapeutics. Trends Pharmacol Sci 2018, 39:59-74.
Wu T, Dai Y: Tumor microenvironment and therapeutic response. Cancer Lett 2017, 387:61-68.
Ge P, Wang W, Li L, Zhang G, Gao Z, Tang Z, Dang X, Wu Y: Profiles of immune cell infiltration and immune-related genes in the tumor microenvironment of colorectal cancer. Biomed Pharmacother 2019, 118:109228.
Bhattacharya S, Andorf S, Gomes L, Dunn P, Schaefer H, Pontius J, Berger P, Desborough V, Smith T, Campbell J, et al: ImmPort: disseminating data to the public for the future of immunology. Immunologic research 2014, 58:234-239.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK: Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun 2019, 10:1523.
Obuchowski NA, Bullen JA: Receiver operating characteristic (ROC) curves: review of methods with applications in diagnostic medicine. Phys Med Biol 2018, 63:07tr01.
Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, Khodadoust MS, Esfahani MS, Luca BA, Steiner D, et al: Determining cell type abundance and expression from bulk tissues with digital cytometry. Nature Biotechnology 2019, 37:773-782.
Zalocusky KA, Kan MJ, Hu Z, Dunn P, Thomson E, Wiser J, Bhattacharya S, Butte AJ: The 10,000 Immunomes Project: Building a Resource for Human Immunology. Cell reports 2018, 25:513-522.e513.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, et al: Inferring tumour purity and stromal and immune cell admixture from expression data. Nature Communications 2013, 4:2612-2612.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA: Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015, 12:453-457.
Liu Y, Lu J, Zhang Z, Zhu L, Dong S, Guo G, Li R, Nan Y, Yu K, Zhong Y, Huang Q: Amlexanox, a selective inhibitor of IKBKE, generates anti-tumoral effects by disrupting the Hippo pathway in human glioblastoma cell lines. Cell Death Dis 2017, 8:e3022.
Ohmori K, Hayashi K, Kaise T, Ohshima E, Kobayashi S, Yamazaki T, Mukouyama A: Pharmacological, pharmacokinetic and clinical properties of olopatadine hydrochloride, a new antiallergic drug. Jpn J Pharmacol 2002, 88:379-397.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A: Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018, 68:394-424.
Court KA, Hatakeyama H, Wu SY, Lingegowda MS, Rodríguez-Aguayo C, López-Berestein G, Ju-Seog L, Rinaldi C, Juan EJ, Sood AK, Torres-Lugo M: HSP70 Inhibition Synergistically Enhances the Effects of Magnetic Fluid Hyperthermia in Ovarian Cancer. Mol Cancer Ther 2017, 16:966-976.
Yang Z, Zhuang L, Szatmary P, Wen L, Sun H, Lu Y, Xu Q, Chen X: Upregulation of heat shock proteins (HSPA12A, HSP90B1, HSPA4, HSPA5 and HSPA6) in tumour tissues is associated with poor outcomes from HBV-related early-stage hepatocellular carcinoma. Int J Med Sci 2015, 12:256-263.
Yao DF, Dong ZZ, Yao M: Specific molecular markers in hepatocellular carcinoma. Hepatobiliary Pancreat Dis Int 2007, 6:241-247.
Wang X, Sun Z, Tian W, Piao C, Xie X, Zang J, Peng S, Yu X, Wang Y: S100A12 is a promising biomarker in papillary thyroid cancer. Sci Rep 2020, 10:1724.
Gunaldi M, Okuturlar Y, Gedikbasi A, Akarsu C, Karabulut M, Kural A: Diagnostic importance of S100A9 and S100A12 in breast cancer. Biomed Pharmacother 2015, 76:52-56.
Cai H, Ye BG, Ao JY, Zhu XD, Zhang YY, Chai ZT, Wang CH, Sun HC: High expression of S100A12 on intratumoral stroma cells indicates poor prognosis following surgical resection of hepatocellular carcinoma. Oncol Lett 2018, 16:5398-5404.
Song GX, Shen YH, Liu YQ, Sun W, Miao LP, Zhou LJ, Liu HL, Yang R, Kong XQ, Cao KJ, et al: Overexpression of FABP3 promotes apoptosis through inducing mitochondrial impairment in embryonic cancer cells. J Cell Biochem 2012, 113:3701-3708.
Tang Z, Shen Q, Xie H, Zhou X, Li J, Feng J, Liu H, Wang W, Zhang S, Ni S: Elevated expression of FABP3 and FABP4 cooperatively correlates with poor prognosis in non-small cell lung cancer (NSCLC). Oncotarget 2016, 7:46253-46262.
Yu LR, Cao Z, Makhoul I, Daniels JR, Klimberg S, Wei JY, Bai JP, Li J, Lathrop JT, Beger RD, Todorova VK: Immune response proteins as predictive biomarkers of doxorubicin-induced cardiotoxicity in breast cancer patients. Exp Biol Med (Maywood) 2018, 243:248-255.
Wang X, Wang SS, Zhou L, Yu L, Zhang LM: A network-pathway based module identification for predicting the prognosis of ovarian cancer patients. J Ovarian Res 2016, 9:73.
Sun W, Shang J, Zhang J, Chen S, Hao M: Correlations of DKK1 with incidence and prognosis of breast cancer. J buon 2019, 24:26-32.
Igbinigie E, Guo F, Jiang SW, Kelley C, Li J: Dkk1 involvement and its potential as a biomarker in pancreatic ductal adenocarcinoma. Clin Chim Acta 2019, 488:226-234.
Stephens JM, Elks CM: Oncostatin M: Potential Implications for Malignancy and Metabolism. Curr Pharm Des 2017, 23:3645-3657.
Tawara K, Scott H, Emathinger J, Wolf C, LaJoie D, Hedeen D, Bond L, Montgomery P, Jorcyk C: HIGH expression of OSM and IL-6 are associated with decreased breast cancer survival: synergistic induction of IL-6 secretion by OSM and IL-1β. Oncotarget 2019, 10:2068-2085.
Feng Q, Wu X, Li F, Ning B, Lu X, Zhang Y, Pan Y, Guan W: miR-27b inhibits gastric cancer metastasis by targeting NR2F2. Protein Cell 2017, 8:114-122.
Zhang C, Han Y, Huang H, Qu L, Shou C: High NR2F2 transcript level is associated with increased survival and its expression inhibits TGF-β-dependent epithelial-mesenchymal transition in breast cancer. Breast Cancer Res Treat 2014, 147:265-281.
Giraldo NA, Sanchez-Salas R, Peske JD, Vano Y, Becht E, Petitprez F, Validire P, Ingels A, Cathelineau X, Fridman WH, Sautès-Fridman C: The clinical role of the TME in solid cancer. Br J Cancer 2019, 120:45-53.
Speiser DE, Ho PC, Verdeil G: Regulatory circuits of T cell function in cancer. Nat Rev Immunol 2016, 16:599-611.
Fridman WH, Pagès F, Sautès-Fridman C, Galon J: The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer 2012, 12:298-306.
Gu-Trantien C, Loi S, Garaud S, Equeter C, Libin M, de Wind A, Ravoet M, Le Buanec H, Sibille C, Manfouo-Foutsop G, et al: CD4⁺ follicular helper T cell infiltration predicts breast cancer survival. J Clin Invest 2013, 123:2873-2892.

Additional file 1: Figure S1. Identification of differentially immune genes in ESCA. a The heatmap and b volcano map of differentially expression genes (n=4102). c Venn diagram was performed to indicate the screening process of differentially immune-related genes. d The heatmap and e Volcano map of 303 differentially immune genes of differentially immune genes.

Additional file 2: Figure S2. Functional enrichment analyses of differentially immune genes. a Gene Ontology (GO)categories and b Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis of differentially immune genes (n=303).

Additional file 3: Table S1. Significant GO terms of differentially immune genes (n=303) in ESCA.

Abbreviations: ESCA, esophageal cancer; GO, Gene Ontology; MF, molecular functions; BP, biological processes; CC, cellular components.

Additional file 4: Table S2. Significant KEGG pathways of differentially immune genes (n=303) in ESCA.

Abbreviations: ESCA, esophageal cancer; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Additional file 5: Figure S3. GO and KEGG pathways of 303 differently immune genes in a network way. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Additional file 6: Figure S4. Assessment the relationship between prognostic signature with infiltrating-immune cells. a-f Systematic association analyses between six immune genes in prognosis signature and several importantly immune cells.

Download PDF

Journal Publication

published 01 Oct, 2020

Read the published version in International Immunopharmacology →

Version 1

posted

You are reading this latest preprint version

Identification the prognostic value of immune gene signature and infiltrating immune cells of esophageal cancer patients

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Materials And Methods

Data acquisition

Identification of differentially immune genes

Construction of the immune-related risk score system

Assessment the prognosis value of immune signature

Exploration the prognosis value of infiltrating immune cells in TME

Statistical analysis

Results

Identification of differentially immune genes in patients of ESCA

Construction of the immune-related prognosis signature in the entire TCGA set

Validation of the predictive performance of immune prognosis signature

Independent prognostic value of six immune-gene signature and clinical parameters

Revealing the relationship between prognostic signature and tumor-infiltrating immune cells in TME

Discussion

Conclusion

Abbreviations

Declarations

References

Supplemental Files

Supplementary Files

Status:

Journal Publication

Version 1