Construction and Evaluation of a Risk Score Model for Autophagy-Related Genes in Esophageal Adenocarcinoma

Background

Esophageal adenocarcinoma (EAC) is one of the most common malignant tumors and the sixth leading cause of cancer deaths. EAC is a pathological subtype of esophageal cancer mainly occurring in the lower esophagus, which is connected with the stomach. EAC development is related to Barrett’s esophagus, gastric reflux, and obesity [1]. At present, the incidence of esophageal adenocarcinoma is increasing in Asia [2]. Owing to the lack of reliable early diagnostic methods and effective interventions, the 5-year survival rate of EAC is less than 20% [3]. EAC malignancy is asymptomatic in the early stage, with symptoms generally appearing in the middle and late stages [4]. Patients with localized EAC can undergo surgical and postoperative therapies, including chemoradiotherapy [5]. However, the treatment of advanced or metastatic EAC provides only palliative and supportive care. Therefore, there is an urgent need to find innovative methods and biomarkers for early detection of EAC and new treatment strategies. So far, several prognostic studies have used various indicators (including DKK3 [6], TXNIP [7], CTO1 [8], and GRB7 [9]) that have been reported as prognostic biomarkers for patients with EAC. However, the median survival of EAC has increased only slightly, and there is still a lack of sensitive gene combination markers to evaluate the prognosis of patients with EAC. Therefore, more studies are needed to explore potential prognostic markers and treatment strategies for EAC.

Autophagy is a crucial biological process that balances cell homeostasis through the degradation of proteins and organelles in cells. Previous studies have validated the role of autophagy in promoting tumor cell survival and suppressing oncogenesis [10–12]; therefore, enhancing and inhibiting autophagy has been suggested as a possible treatment strategy [13–16]. Other studies have demonstrated that autophagy activation shows a marked correlation with tumor dormancy, chemoresistance, and stem cell survival [17]. Exploring the relationship between EAC and autophagy could further elucidate the mechanisms of pathogenesis and improve prognosis and treatment options.

Large-scale expression data have scarcely assessed the role of autophagy in EAC prognosis. This study constructed a bioinformatic prognostic model of EAC based on the RNA-seq expression profile of EAC tissues in The Cancer Genome Atlas (TCGA) database and further determined predictive biomarkers related to autophagy.

Material and Methods

PATIENT SAMPLES AND GENE EXTRACTION:

EAC-related RNA-seq data and clinical data were downloaded from the TCGA database (TCGA, https://portal.gdc.cancer.gov/). EAC clinical data files included survival status, sex, age, and TNM staging. Perl language (http://www.perl.org/) was used to combine RNA-seq data into a single matrix file and normalize the matrix data. In the process, if multiple probes detected the expression level of the same gene, the average expression level of the gene was taken as the expression value. In the data processing of patients’ clinical data, the clinical information of patients whose survival status was unknown and whose survival time was blank was deleted. Clinical data included for statistical analysis included patient number, survival time, survival status, age, sex, classification, and TNM stage. The local ethics committee approved all open-access datasets and related research, and our further data mining of these datasets complied with the principles of the Declaration of Helsinki.

A query in the human autophagy database (HADb, http://www.autophagy.lu/index.html) showed that, to date, there are 232 autophagy-related genes (ARGs). EAC data and ARGs were obtained by screening RNA-seq data downloaded from the TCGA database.

UNIVARIATE COX AND LASSO REGRESSION ANALYSIS:

Using combined patient survival information, the survival package in R software was used to conduct a univariate Cox regression analysis of genes with autophagy-related differences. The hazard ratio (HR) and P value of each gene and the survival of patients with EAC were calculated, and the candidate genes significantly related to the prognosis of patients with EAC were screened out according to the criterion of P<0.05.

The screened genes were further analyzed by lasso regression analysis using the glmnet package of R software [18,19], and the regression coefficient was determined by cross-validation to reduce the influence of collinearity among genes [20], prevent the risk model variables from overfitting, and conduct other screens affecting the prognosis of EAC.

MULTIVARIATE COX REGRESSION ANALYSIS AND PROGNOSTIC MODEL CONSTRUCTION:

Using survival R software packages, the lasso method was used to select genes obtained by multivariable Cox regression analysis. According to the Akaike information criterion (AIC), the lowest AIC value yields the smallest model, which is the optimal model [21] containing the critical genes. To calculate the critical genes multifactor regression coefficient, a prognostic risk factor score equation was built:

where exp is the expression level of each gene in the sample, and β is the corresponding gene’s multifactor regression coefficient. The prognostic risk coefficient of each gene was calculated according to the formula. According to the risk value’s median value, the HR value and 95% confidence interval (CI) of the critical genes in the model were calculated. Based on the median risk score, patients with EAC were divided into low- and high-risk groups. The survival analyses were recognized as key outcomes, and Kaplan-Meier (K-M) analysis was used to analyze the overall survival among the various score ranges. Time-dependent receiver operating characteristic (ROC) curves were plotted to calculate the risk model’s predictive ability. R software (version 3.6.1) was used to make bioinformatic analyses. P<0.05 was regarded as statistically significant.

PROGNOSTIC MODEL EVALUATION AND ONCOMINE ANALYSIS:

The ROC curve of the model was drawn using the survival ROC package in R software. We determined the area under the curve (AUC) values, and the C-index was calculated to evaluate the model’s sensitivity and prediction ability. The higher the AUC and C-index, the better the model’s prognostic accuracy. Through K-M survival analysis of the critical genes in the prognostic model, the value and significance of a single gene in EAC prognostic prediction were explored.

The differentially expressed genes were analyzed using the Oncomine database (www.oncomine.org) [22]. The prognostic model results were verified, and the search conditions were set as high-expression ATG5 gene and low-expression BECN1 gene in the survival analysis of the essential prognostic genes. For differential analysis, the analysis parameters were P<0.05, fold change=all, gene rank=all, cancer type=esophageal adenocarcinoma, and the first 4 results with the lowest P values were drawn.

GENE SET ENRICHMENT ANALYSIS:

To further explore the biological functions of EAC disease-related genes, gene set enrichment analysis (GSEA) was used in this study [23] to assess the phenotypic correlation gene in the distribution trend list, sort and analyze critical genes in EAC prognosis, and identify the biological function of a participating signaling pathway. From the Molecular Signature database (MsigDB), the Gene Ontology (GO) reference gene set “c5. all. v7.1. symbols. gmt” and Kyoto Encyclopedia of Genes and Genomes (KEGG) reference gene set “c2. cp. kegg. v7.1. symbols. gmt” were used. Enrichment analysis was carried out using the weighted enrichment method. The random combination times were set as 1000, and the items with the highest enrichment scores and statistical significance (P<0.05) were visualized for analysis.

Results

SCREENING FOR DIFFERENTIAL GENES RELATED TO AUTOPHAGY:

The RNA-seq expression profile and clinical data of EAC were screened by the TCGA database, including 78 EAC samples and 9 normal tissue samples. The clinical baseline data of EAC patients are summarized in Table 1. The HADb database was used to screen 206 ARGs.

UNIVARIATE COX AND LASSO REGRESSION ANALYSIS:

A univariate Cox regression analysis was performed on ARGs to calculate the HR and P value of each gene for survival in EAC patients; 14 genes with a significant difference (P<0.05) were screened. Significant association defines the gene as a candidate gene, and the genes were mapped on a single factor Cox regression analysis diagram (Figure 1).

These genes were further screened using lasso regression (P<0.05) (Figure 1), and a model was established using cross-validation (Figure 2) to screen 13 significantly correlated candidate genes.

MULTIVARIATE COX ANALYSIS:

Thirteen significantly correlated candidate genes were screened by multivariate Cox analysis and the optimal model was selected according to the lowest AIC value; 8 genes were defined as key in the model and were used to draw the gene multivariable Cox analysis diagram (Figure 3).

The prognostic risk score was calculated according to the corresponding regression coefficient, where the risk score=(0.460965×TP73exp)+(−0.4529×DAPK1exp)+ (−1.1353×BECN1exp)+(1.75806×ATG12exp)+ (0.680104×SIRT1exp)+(−0.92933×CAPN1exp)+ (1.647224×VAMP7exp)+(1.638798×ATG5exp).

The prognostic risk index of all samples was calculated and divided into cutoff points by the median value dividing the risk indexes into a high-risk group and a low-risk group for subsequent analysis.

DRAWING OF SURVIVAL CURVE, RISK CURVE, AND ROC CURVE:

The survival score was obtained based on a prognostic risk assessment model derived from the multivariate Cox analysis. The median of the prognostic risk index of the sample was 1.14E+00. Patients with a risk index >1.14E+00 were defined as a high-risk group and patients with a risk index <1.14E+00 were defined as a low-risk group.

The risk curve, survival point map, and risk value of the relationship between the expression of critical genes were drawn by analyzing the survival difference of patients in the high- and low-risk assessment groups (Figure 4A). The results showed that with the increase of risk score, the number of deaths in the high-risk group was higher than that in the low-risk group, indicating that the survival rate of patients in the high-risk group was low. A K-M survival curve (Figure 4B) was drawn for the high-risk and low-risk groups. From the survival data, the 3-year survival rate was 13.1% (95% CI, 4.15–4.16%) in the high-risk group and 70.2% (95% CI, 53.3–92.6%) in the low-risk group. The difference in survival rate between the 2 groups increased with the passing of time. The difference between the 2 groups was statistically significant (P<5.978E-06). At the same time, the ROC curve model was drawn (Figure 4C), where the overall survival for AUC was 0.849, and the C-index was 0.808, indicating that the model was accurate in the prognosis of patients with EAC and had good clinical value.

SURVIVAL ANALYSIS OF ESSENTIAL EAC PROGNOSTIC GENES:

Further single-gene K-M survival analysis was performed on the 8 critical genes in the model obtained by multivariate Cox regression analysis. The results showed that 4 of the genes were statistically significant (P<0.05) and related to prognosis, namely, VAMP7, ATG5, BECN1, and DAPK1. The higher the expression level of VAMP7 and ATG5, the worse the prognosis, and the higher the expression level of BECN1 and DAPK1, the better the prognosis (Figure 5). These results indicated that VAMP7 and ATG5 may be independent adverse factors for EAC prognosis, while BECN1 and DAPK1 may be independent protective factors for EAC prognosis.

DIFFERENTIAL EXPRESSION ANALYSIS OF ATG5 AND BECN1 GENES IN THE ONCOMINE DATABASE:

The Oncomine database was used to analyze the differential expression of ATG5, an independent adverse prognostic factor of EAC, and BECN1, an independent protective factor of EAC. The first 4 differential analyses results of the 2 genes with statistical significance (P<0.05) were selected for statistical analysis. Among all differentially expressed genes, the median value of ATG5 was 580.5 (P=0.008), indicating that ATG5 was significantly highly expressed in EAC tissue (Figure 6A). The BECN1 gene ranked 1739.5 among all differentially expressed genes (P=0.007), indicating that BECN1 had significantly low expression in EAC tissue (Figure 6B).

GSEA OF ESSENTIAL PROGNOSTIC GENES:

GO functional enrichment analysis of prognostic EAC genes by GSEA was carried out with “c5. all. v7.1. symbols. gmt” as the reference gene set, and the critical prognosis genes in GO were enriched in 5426 expression samples, 3949 of which had high expression and 1477 of which had low expression. There were 168 statistically significant genes in the high-expression samples (P<0.05) and 28 statistically significant genes in the low-expression samples (P<0.05). The top 5 items with highest enrichment scores in the high- and low-expression sample (Table 2) were visualized (Figure 7).

Using “c2. cp. kegg. v7.1. symbols. gmt” as the reference gene set, KEGG enrichment analysis results showed the prognosis of critical gene enrichment of the KEGG pathway in 178 genes, 94 of which had high expression and 84 of which had low expression. Among them, 3 were statistically significant in the high-expression samples (P<0.05), and 6 were statistically significant in the low-expression samples (P<0.05) (Table 3). The above 9 items were visualized in high and low expression (Figure 8).

Discussion

Although esophagectomy, chemotherapy, and radiotherapy can improve EAC prognosis, not all EAC patients are helped by these therapies. The current application of molecularly targeted drugs in the treatment of EAC is not satisfactory [24]. The mechanism of autophagy in tumors has long been proposed as a new therapeutic approach. Several studies have supported an association of distinct expression levels between autophagic genes and EAC, showing that ARG regulation may be associated with EAC occurrence and development [25]. Therefore, we searched for EAC autophagy genes as molecular markers to help design an autophagy-targeted treatment strategy to improve patient prognosis with esophageal cancer.

This study summarized the RNA-seq data and clinical information on EAC in the TCGA database. The HADb database screened 206 ARGs, univariate Cox regression analysis was used to obtain 14 prognostic ARGs to reduce the collinearity between genes, 13 ARGs were obtained by Lasso regression analysis, and the HR values and corresponding P values of autophagy-related prognostic genes were further calculated by multivariate Cox regression analysis. Eight prognostic genes related to EAC and autophagy were screened and constructed. Those 8 genes were defined as critical genes in the risk score model of EAC and ARGs, and the prognostic risk model was constructed. The median of the prognostic risk index of the model was 1.14E+00. According to the prognostic risk index’s median value, patients with EAC were divided into high-risk and low-risk assessment groups. The model’s K-M risk curve was drawn, the ROC curve of the model was drawn, and the AUC was calculated as 0.849. The C-index was 0.808, and the model’s AUC and C-index showed that the model was accurate for the prognosis of EAC in patients, with good clinical practice value to determine the diagnostic and prognostic value of essential prognostic ARGs in EAC. K-M survival analysis of the essential prognostic genes in the multivariate Cox regression model was performed. The obtained 4 genes, DAPK1, BECN1, ATG5 and VAMP7, showed statistical significance (P<0.05) and were related to prognosis. The higher the expression level of VAMP7 and ATG5, the worse the prognosis, and the higher the expression level of BECN1 and DAPK1, the better the prognosis. ATG5 and BECN1 were selected for Oncomine differential expression analysis, showing that ATG5 was significantly highly expressed in EAC tissue. The expression of BECN1 was significantly low in EAC tissue, which verified the survival analysis results of the critical gene by K-M. In GSEA with “c5” as the reference gene set, it was found that the high and low expression of the gene set was mainly concentrated in biological functions such as cell generation and regulation of the cell cycle. In GSEA with “c2” as the reference gene set, it was found that the high and low expression of the gene set was mainly concentrated in metabolic pathways such as maintaining cell activity. Four genes associated with prognosis and with statistical significance (P<0.05) were analyzed.

DAPK1 regulates multiple cellular processes and plays an essential role in the pro-apoptotic, apoptotic, and autophagy pathways. DAPK1 participates in autophagy by phosphorylating downstream Beclin-1 and enhances autophagy via interaction with Beclin-1 [26]. Some studies have found that Beclin-1 gene mutation or gene knockout cells could not form an effective autophagosome to clear damaged organelles, and mitochondria could not repair DNA and protein damage, leading to aggregation of harmful substances, inflammatory factors, and cytokines, thus inducing tumor formation [27]. Other studies showed that Beclin-1 expression is significantly decreased in the tumor tissues of digestive tract tumors and the corresponding normal tissue cells, indicating that Beclin-1-mediated autophagy plays a role in the inhibition of tumor development [28]. Recent research shows that the loss of Beclin-1 expression occurs more frequently in patients with EAC than in controls (49.0% vs. 4.8%) [29]. Additionally, Beclin-1 expression level is negatively correlated with EAC histologic grade and stage (P<0.005). Beclin-1 could cause long-term survival, which implies that Beclin-1 is a tumor suppressor gene that may act as a prognostic biomarker. The above research results are consistent with the results of the present study. The autophagic protein ATG5 plays an essential role in the formation of the autophagic vacuole [30], which could combine with Atg12 and Atg16 to form an autophagic complex [14]. Studies have found that increased ATG5 expression was observed in esophageal cancer, compared to normal tissue, and the high ATG5 expression tumor group showed a poor prognosis, including overall survival and progression-free survival [31]. In addition, ATG5 inhibition contributed to the treatment of patients with esophageal carcinoma [32]. This result was consistent with the results of our present study. Another study found that ATG5 protein expression in several gastric cancer cell lines was higher than that in regular gastric mucosa epithelial cell lines; the high expression of ATG5 in gastric cancer promoted the occurrence and development of tumors [33]. The above studies show that ATG5 is a high-risk gene, which indicates worse prognosis when highly expressed. VAMP7 is a SNARE protein, mediating specific membrane fusion through intracellular transport. It has been reported that VAMP7 can regulate the formation of autophagy [34]. Some studies have found that VAMP7 connects with the tumor cell membrane through protease secretion vesicles, and that the protease accumulated at the end of pseudopodia erodes the basement membrane surrounding breast cancer tissue [35]. VAMP7 inactivation could significantly reduce the ability of breast cancer cells to degrade the extracellular matrix, which may be the molecular mechanism of potential tumor cell metastasis. These prognostic ARGs may be useful for cancer detection and might be a valid strategy to increase patient survival.

Conclusions

ARGs affecting prognosis were screened and combined into a risk scoring model. Multiple bioinformatic methods confirmed their essential role in EAC prognosis. The relationship between risk score and survival time of patients was verified. The expression of these 4 genes, DAPK1, BECN1, ATG5 and VAMP7, can be used to calculate the risk score for any patient with EAC to evaluate the patients by clinical data analysis. The model could also be used as an independent prognostic factor, providing adequate prognostic information and clinical application value for patients with EAC.