Identification of DEGs
After data preprocessing, the RNA SEQ data of 159 EC specimens and 11 adjacent specimens were included. The volcano map (Fig. 1A) shows that 1220 DEGs are differentially expressed between EC and adjacent samples, in which the red dot represents significantly up-regulated genes, the green dot represents significantly down-regulated genes, and the black dot represents no difference genes.
Differential gene analysis and prognosis analysis of esophageal cancer. (A) 1220 DEGs differentially expressed between the EC and adjacent tissues. Red dots represent differentially expressed up-regulated genes, green dots represent differentially expressed down-regulated genes, and black dots represent no significant difference in gene expression. (B) Kaplan–Meier (KM) prognostic analysis of DEGs. 41 genes of 1220 genes in the DEGs were associated with the prognosis of patients with esophageal cancer.
Screening of prognostic genes in EC
Kaplan-Meier (KM) survival analysis was used to analyze the proteins related to the prognosis of the EC. The samples were divided into high expression or low expression (relative to median expression). The survival time and status of low expression group and high expression group were compared. Among 1220 DEGs, 41 proteins were associated with the prognosis of the EC (Fig. 1B). Among them, 22 genes were up-regulated and 19 genes were down regulated (Fig. 2).
Up and down regulation of prognostically related DEGs in esophageal cancer. The red bar represents the up regulation of gene expression, the blue bar represents the down regulation of gene expression, and the length of the bar represents the value of |logFC|. The prognostic p value was shown on the right side of the gene.Up-regulated:APOC1,TMEM270,VWDE,APOA2,SIX3,CSF2,TAS2R38,YBX2,METTL27,PAEP,HIST1H2AJ,SOST,SHISA2,MMP12,CT45A1,HS6ST2,HIST1H2BI,CLDN3,UBE2C,POU6F2. Down-regulated:GPER1,KCTD8,NEXMIF,RYR2,SNAP91,COX6A2,CPA2,FAM189A2,SULT2A1,RBFOX3,CLCNKB,KCNK2,NSG2,BHLHA15,FABP3,KCNG4,CAPZA3,ERO1B,GRIA3.
Functional enrichment analysis
GO and KEGG pathway enrichment analysis was performed on 41 prognostic related DEGs. The results show that the GO analysis shows that the biological process (Fig. 3A) includes phospholipid flux, regulation of cholesterol estimation, glycolipid catabolic process, high-density lipoprotein particle remodelling, steroid estimation, sterol estimation, cholesterol estimation, crathrin coat assembly, negative regulation of lipase activity, Phosphoridylcholine metallic process. In cellular component (Fig. 3B), it includes ion channel complex, chylomicron, transmembrane transporter complex, transporter complex, very low density lipoprotein particle, triglyceride rich plasma lipoprotein particle, cation channel complex, high density lipoprotein particle, plasma lipoprotein particle and lipoprotein particle. The molecular function (Fig. 3C) includes gated channel activity, lipase inhibitor activity, ion channel activity, phosphotylcholine binding, quaternary ammonium group binding, channel activity, passive transmembrane transporter activity, fatty acid binding, location channel activity and sulfotransferase activity. Through GO analysis, it is found that DEGs play a role in lipid metabolism, channel complex composition and channel activity, and gene differential expression may lead to the disorder of the above biological functions of cells. KEGG pathway is mainly concentrated in PPAR signaling pathway, cardiac muscle contract and pancreatic secret (Fig. 3 D). Among them, previous studies believe that peroxisome promoter activated receptor (PPAR) signaling pathway has PPAR signal imbalance in a variety of cancers and focuses on a variety of metabolic pathways14.
GO and KEGG enrichment analysis of 41 DEGs related to prognosis in EC. (A) Biological Process (BP). (B) Cellular Component(CC). (C) Molecular Function (MF). (D) Kyoto Encyclopedia of Genes and Genomes (KEGG). (According to the p value from small to large, the color of the circle increases from red to blue. The size of the circle indicates the count, that is, the number of genes.)
Construction and evaluation of prognostic model
The clinical follow-up data of 159 patients with EC were combined with 1220 DEGs. Univariate Cox analysis showed that 42 DEGs affected the prognosis of EC. Venn diagram shows 15 common prognostic related genes obtained by Kaplan-Meier (KM) and univariate Cox methods (Fig. 4). Then, multivariate Cox regression analysis was performed on these 15 candidate genes, and 8 DEGs (APOA2, COX6A2, CLCNKB, BHLHA15, HIST1H1E, FABP3, UBE2C and ERO1B) were significantly correlated with OS in patients with EC (Fig. 5). Among them, the HR values of APOA2, HIST1H1E, FABP3, UBE2C and ERO1B were greater than 1, which were potential risk factors. The HR values of COX6A2, CLCNKB and BHLHA15 were less than 1, which were potential protective factors. The difference of prognosis of eight genes was calculated by survival curve (Fig. 6). After integrating 8 genes and weighting their multivariable Cox regression coefficients, the risk score formula was obtained: (0.13592497×Expression of APOA2) + (-0.267351021×Expression of COX6A2) + (-0.26668478×Expression of CLCNKB) + (-0.265751714×Expression of BHLHA15) + (0.30430453×Expression of HIST1H1E) + (0.4497437×Expression of FABP3) + (0.336614943×Expression of UBE2C) + (0.264739114×Expression of ERO1B). According to the risk score formula, 159 patients with EC were given a risk value and divided into high-risk group and low-risk group with the median risk value as the cut-off value. The grouping results were visualized by prognostic feature distribution map (Fig. 7A), 8 DEGs expression profile heat map (Fig. 7B), patient survival map (Fig. 7C), ROC curve (Fig. 7D) and survival curve (Fig. 7E). KM survival curve showed that the survival rate of high-risk group was significantly lower than that of low-risk group (P = 8.124e-07).
Venn diagrams for screening prognostic genes of the EC by Kaplan–Meier (KM) and univariate Cox analysis. A total of 15 common prognostic related genes (GPER1,KCTD8,TMEM270,APOA2,COX6A2,CPA2,FAM189A2,CLCNKB,BHLHA15,HIST1H1E,FABP3,HIST1H2BI,CLDN3,UBE2C and ERO1B) were found.
Forest maps of multivariate Cox regression analysis results. HR (Hazard Ratio) represents the risk coefficient of high expression group relative to low expression group. If HR > 1, the gene is a risk factor; If HR < 1, the gene is a protective factor; 95% Cl represents HR confidence interval. * p < 0.05, ** p < 0.01, *** p < 0.001.
Kaplan–Meier (KM) survival curves of the EC patients with low or high individual risk score of eight genes. (A) Survival curve of APOA2. (B) Survival curve of COX6A2. (C) Survival curve of CLCNKB. (D) Survival curve of BHLHA15. (E) Survival curve of HIST1H1E. (F) Survival curve of FABP3. (G) Survival curve of UBE2C. (H) Survival curve of ERO1B. The red line represents the high-risk group and the blue line represents the low-risk group.
Performance evaluations of prognostic risk scoring model. (A) The patient’s risk score, red indicates high risk and green indicates low risk. (B) Expression heat map of eight genes in high-risk group and low-risk group (Blue represents the high-risk group and red represents the low-risk group). (C) Distribution of survival status of patients in high-risk group and low-risk group (red indicates death and green indicates survival). (D) ROC curve of comprehensive risk scores of eight genes. AUC (area under curve) indicates the area below the ROC curve. The value is between 0 and 1. The higher the value, the better the prediction effect of the model. (E) KM survival curve of patients with the EC with low or high comprehensive risk score of eight genes.
Evaluation of prognostic model as an independent prognostic factor
The clinical characteristics of different individuals may affect their prognosis. Therefore, the calculated risk score and other clinical characteristics (age, gender, T, N, M and tumor stage) were included in univariate and multivariate Cox regression analysis. In univariate analysis (Fig. 8A), clinical characteristics (gender, age, T, N, M and tumor stage), APOA2, COX6A2, HIST1H1E, UBE2C, ERO1B and eight gene comprehensive risk scores were prognostic risk factors for patients with the EC (P < 0.05). In multivariate analysis (Fig. 8B-J), individual prognostic risk scores of APOA2 (P = 0.03), COX6A2 (P = 0.012), BHLHA15 (P = 0.022), HIST1H1E (P = 0.019), FABP3 (P = 0.009) and UBE2C (P = 0.017) were significantly correlated with prognosis. At the same time, the eight gene comprehensive risk score showed a stronger prognostic correlation (P < 0.001), indicating that the eight gene comprehensive risk score can be used as an independent predictor. ROC curve analysis is used to evaluate the prediction efficiency (Fig. 8K-S). The comprehensive risk scores of the eight genes under the 1-year, 3-year and 5-year curves were 0.718, 0.862 and 0.95 respectively, which was a better predictor than other characteristic factors.
Evaluation of eight gene models as independent predictors. (A) ROC curve with comprehensive risk scores as a predictor. (B–J) Multivariate analysis of eight gene individuals and comprehensive risk scores involving patient characteristics. (K–S) ROC curve. The comprehensive risk scores of the eight genes under the 1-year, 3-year and 5-year curves were 0.718, 0.862 and 0.95.
Regulatory network of risk gene-miRNA and drug sensitivity analysis
This study further collected and analyzed the miRNA expression profile data of the EC from the TCGA database. Among them, 98 up-regulated and 64 down regulated differentially expressed microRNAs were included (Fig. 9A). MirDIP analysis of the regulatory relationship between differential microRNAs and 8 prognostic mRNAs showed that 12 differentially expressed microRNAs had 13 potential regulatory relationships with 3 prognostic mRNAs (Fig. 9B). The results of drug sensitivity of risk genes showed that the high expression of HIST1H1E made tumor cells resistant to trametinib, selumetinib, RDEA119, Docetaxel and 17-AAG. The high expression of UBE2C makes tumor cells resistant to masitinib. The low expression of ERO1B makes the EC more sensitive to FK866 (Fig. 10).
Construction of gene miRNA regulatory network. (A) Volcano map of differential expression miRNA screening. Green dots represent down regulated miRNAs and red dots represent up regulated miRNAs. (B) Regulatory network between risk genes and differentially expressed miRNAs. Red nodes represent up-regulated miRNAs and blue nodes represent down-regulated miRNAs.
Potential drugs screened from GSCA database. Red nodes represent positive correlation and blue nodes represent negative correlation. The larger the node, the stronger the correlation. Positive correlation means that the high expression of the gene is resistant to drugs, and vice versa.
Expression of UBE2C in tissues and cells
In immunohistochemical analysis, UBE2C was strongly expressed in esophageal cancer tissues (Fig. 11A) and negatively expressed in adjacent tissues (Fig. 11B). The expression was relatively strong in esophageal cancer cell lines (kyse150, TE-1 and Eca109) (Fig. 11C-E). The expression was negative in normal esophageal cell line (HEEC) (Fig. 11F). 150 pairs of cancer and adjacent tissues were concentrated and immunohistochemical expression chips were constructed (Fig. 11G). Subsequent statistical analysis found that the expression of UBE2C in esophageal cancer tissues was higher than that in adjacent tissues (Fig. 11H). ROC curve results showed that UBE2C had a good differential diagnosis ability for esophageal cancer (Fig. 11I). Correlation analysis between UBE2C expression and clinical factors showed that UBE2C expression was associated with tumor metastasis, and the positive rate of metastasis was high in patients with high expression (Table 1).
Immunohistochemical staining of UBE2C in the tissues and cell lines of the EC. (A) Expression of UBE2C in the EC. (B) Expression of UBE2C in normal esophageal tissues. (C) Expression of UBE2C in esophageal cancer cell line kyse150. (D) Expression of UBE2C in esophageal cancer cell line Eca109. (E) Expression of UBE2C in esophageal cancer cell line TE-1. (F) Expression of UBE2C in HEEC (normal esophageal cell line) . (G) Expression of UBE2C in the tissue microarray assay of 150 patients with esophageal cancer. (H) Differential expression analysis of UBE2C in cancer and adjacent tissues in the tissue microarray assay. Paired sample t-test was used to compare cancer and adjacent samples. ***P < 0.001. (I) Subject operating characteristic (ROC) curve analysis and area under curve (AUC) statistics were used to evaluate the ability of UBE2C to distinguish esophageal cancer from adjacent normal tissues.(AUC:0.879).
Correlation between UBE2C and clinical indexes and diagnostic efficacy
According to the results of database analysis, UBE2C was positively correlated with the expression of tumor markers (BRCA1, KI67 and TP53) in esophageal cancer (Fig. 12A,D,G). According to the analysis of clinical data, the expression of tumor markers (BRCA1, KI67 and TP53) is strong in cancer tissues (Fig. 12B,E,H). At the same time, UBE2C was positively correlated with the expression of clinical markers (BRCA1, KI67 and TP53) (Fig. 12C,F,I). ROC curve analysis shows that the areas under the curve of BRCA1, KI67 and TP53 are 0.927, 0.940 and 0.902 respectively (Fig. 12J,K,L). However, UBE2C combined with clinical markers (BRCA1, KI67 and TP53) calculated the highest area under the ROC curve, which was 0.996 (Fig. 12M).
Correlation analysis between UBE2C and clinical markers and evaluation of combined diagnostic effect. (A, D, G) The correlation between UBE2C and clinical markers (BRCA1, KI67 and TP53) in esophageal cancer was analyzed by GEPIA database. R represents the correlation coefficient. (B, E, H) Expression of clinical markers (BRCA1, KI67 and TP53) in esophageal cancer. The brown part represents the target protein. (C, F, I) The correlation between UBE2C and clinical markers (BRCA1, KI67 and TP53) was analyzed according to the results of clinical immunohistochemistry. (J-L) ROC curve was used to analyze the diagnostic efficacy of clinical markers (BRCA1, KI67 and TP53) in esophageal cancer. (M) Efficacy of UBE2C combined with clinical markers (BRCA1, KI67 and TP53) in the diagnosis of esophageal cancer.
