Identification and validation of an autophagy-related gene signature for predicting prognosis in patients with esophageal squamous cell carcinoma

RNA expression profiles and corresponding clinical information in tumor and normal tissues of 179 ESCC patients are publicly available (GSE53625↗)¹⁰. The GSE53625↗ dataset contained normalized gene expression data, but the batch effect was still present, so the “limma” package was used to remove the batch effect (Supplementary Fig. 1)¹¹. The maximum expression was considered to be the expression of genes with several probes. 232 and 572 autophagy-related genes were searched from Human Autophagy Database (HADb, http://www.autophagy.lu/↗) and Molecular Signatures Database (version 7.2, http://www.gsea-msigdb.org/gsea/index.jsp↗), respectively. Finally, the expression levels of 585 ARGs were obtained from GSE53625↗ dataset (Supplementary Table 1).

By comparing paired tumor to normal tissues from 179 ESCC patients, differentially expressed autophagy-related genes (DE-ARGs) were identified by absolute value of log2 fold change (log2FC) > 1 and false discovery rate (FDR) < 0.05.

In order to explore the potential mechanism of DE-ARGs affecting ESCC, the “clusterProfiler” package was used to analyze the function and pathway enrichment of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG, www.kegg.jp/kegg/kegg1.html↗)^12–17.

All ARGs were analyzed by univariate Cox regression analysis, and genes with P < 0.05 were defined as prognosis-related genes. Then, multivariate Cox analysis was performed based on the Akaike Information Criterion (AIC) to establish the optimal prognostic risk model. Furthermore, the ARGs-related prognostic model was established by multiplying multivariate Cox regression coefficients by the expression level of each variable.

The risk score of ESCC patients was calculated by using the prognostic signature. According to the median risk score, all ESCC patients were separated into high-risk group and low-risk group. The difference of overall survival between the two groups was analyzed by Kaplan–Meier survival curve. Receiver operating characteristic (ROC) curve and area under ROC (AUC) were used to evaluate the sensitivity and specificity of the prognostic signature. The Brier score and calibration curves were used to assess the calibration of the prognostic model.

In order to verify the effectiveness of the prognostic model, the bootstrap method based on 1000 resampling was used to obtain the testing set^18,19. The training set was the original dataset. In the testing set, AUC and Brier score for the prognostic model were calculated using the R package of the “riskRegression”. The relationship between the risk score and clinical factors was also analyzed to identify the validity of the risk signature, and the survival prediction ability of prognostic factors was further compared. By comparing the clinical traits, univariate and multivariate Cox regression analyses were used to confirm the independence of the prognostic model. At the same time, a nomogram was constructed using the Cox regression coefficients with the R package “rms”, and its calibration curves were drawn.

GSEA, using the “clusterProfiler” package, was employed for assessing the possible mechanisms between high-risk and low-risk groups of ESCC patients based on the KEGG (v7.2) and Hallmark (v7.2) gene sets collections.

The CIBERSORT algorithm was employed to calculate the relative proportion of 22 tumor-Infiltrating Immune Cells (TICs) in 179 ESCC samples^20,21. After quality screening (P < 0.05), 107 cases of ESCC could be for subsequent analysis. The Spearman coefficient examined the correlations between risk score and 22 TICs proportion. Meanwhile, Spearman correlation was also applied to detect the relationship between the risk score and immunotherapy-related targets, such as programmed cell death 1(PDCD1, also known as PD-1) and programmed cell death 1 ligand 1(PD-L1, also known as CD274).

Ferroptosis is a newly introduced form of programmed cell death discovered in recent years. Accumulating studies have revealed ferroptosis could inhibit tumor growth or promote tumor proliferation in tumor development, and there is crosstalk with autophagy at the molecular level^22,23. Then, Spearman correlation examined the correlation between the risk score and ferroptosis -associated genes, which were identified from the FerrDb website (http://www.zhounan.org/ferrdb/index.html↗).

To minimize the bias, the expression of ARGs in the prognostic signature was detected using the following databases: The Cancer Genome Atlas (TCGA) (http://cancergenome.nih.gov/↗), Oncomine²⁴ and Cancer Cell Line Encyclopedia²⁵. RNA-seq data files from TCGA ESCC samples were analyzed as previously described²⁶. Kaplan–Meier plotter database (http://www.kmplot.com/↗) was utilized to validate the prognosis of critical ARGs in the prognostic signature²⁷.

Human ESCC cell lines (Eca-109, KYSE-150 and TE-1) were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Human normal esophageal cell line (Het-1A) was purchased from ATCC (American Type Culture Collection, Manassas, USA). Cell culture, total RNA extraction and qRT-PCR analysis were carried out as described previously²⁸. The relative expression of each gene in prognostic signature was calculated after normalization to GAPDH. Primer sequences are listed in Supplementary Table 2.

Paraffin-embedded ESCC tissues (n = 42) were obtained from pathologically proved ESCC patients who received esophagectomy in the First Affiliated Hospital of Xi'an Jiaotong University from January 2010 to January 2011. In addition, 10 cases of paraffin-embedded normal esophageal mucosa were used as controls. The experiments, including any relevant details, were approved by the Ethics Committee of the First Affiliated Hospital of Xi’an Jiaotong University. The study was performed in accordance with relevant guidelines and regulations. Verbal informed consent was obtained from the patient(s) for their anonymized clinical information to be published in this article. The detailed clinical data of ESCC patients are supplied in Supplementary Table 3. IHC staining analysis was conducted as described previously²⁹. The primary antibody used was anti-cortactin (CTTN) (ab81812, Abcam, 1:300).

Statistical analyses were performed in R software (Version 4.0.2; https://www.R-project.org↗)³⁰. Statistical significance was set at P-value < 0.05. The Kaplan–Meier method was used for survival analysis, and the log-rank test was used to compare the results. Wilcoxon test was used for two independent samples or the paired data. One-way analysis of variance (ANOVA) or Kruskal–Wallis test was used to test the difference among multiple groups. P < 0.05 and Benjamini–Hochberg adjusted P < 0.05 were considered as the cut-off criterion for enrichment analysis.

The work was approved by the Ethics Committee of The First Affiliated Hospital of Xi'an Jiaotong University.

The calibration curves of prognosis signature for the prediction of 1-, 3- and 5-year OS rates demonstrated good agreement in the training set, with Brier scores of 0.15, 0.197, and 0.175, respectively (Fig. 4D). Through bootstrapping verification, the AUCs and brier scores of predicted 1-, 3- and 5-year OS rates were 0.714, 0.714, 0.751, 0.164, 0.222 and 0.203, respectively, indicating that the model had good discrimination and calibration (Fig. 4E).

The analysis of the risk score and clinicopathological features showed that the risk score was associated with survival outcomes (P < 0.001), stage (P = 0.033) and N stage (P = 0.002) (Fig. 4F-H, Supplementary Fig. 2). Furthermore, Fig. 4I showed the 5-year AUC value of the prognostic signature was significantly higher than that of age (AUC = 0.561), gender (AUC = 0.459), grade (AUC = 0.581), stage (AUC = 0.622) and T stage (AUC = 0.518) and N stage (AUC = 0.646).

In order to verify the accuracy of bioinformatics analysis, we used ESCC cell line and ESCC tissue to detect the expression level of ARGs in prognostic signature. The qRT-PCR results were consistent with the above bioinformatics analysis. Compared with human normal esophageal epithelial cells (Het-1A), CTTN was increased and CFTR was down-regulated in ESCC cell lines (Fig. 9G, H). At the cellular level, the expression patterns of the other six genes were similar to those of previous bioinformatics analysis (Supplementary Fig. 6). To further study the CTTN protein expression, we performed immunohistochemistry staining of 42 ESCC tumor samples and 10 normal esophageal samples. The results indicate the expression of CTTN in tumor tissues was significantly higher than that of normal esophageal tissues (Fig. 9I, Supplementary Table 11). Moreover, Kaplan–Meier survival analysis showed that OS and progression-free survival were shorter in the group with high CTTN than in the group with low- negative CTTN (Fig. 9J, K).

Supplementary Information.