Identification of a nomogram based on long non-coding RNA to improve prognosis prediction of esophageal squamous cell carcinoma

DELs identification

All 502 lincRNAs was identified as differentially expressed lincRNAs (DELs) between tumor tissues and adjacent normal tissues, including 223 up-regulated and 279 down-regulated genes (Figure 1A and 1B).

Figure 1. Identification of differentially expressed lincRNAs (DELs) by using “edgeR” and “DEseq” R package. (A) The heatmap of the DELs in ESCC when compared with normal tissue. (B) Volcano plot shown the expression change in ESCC when compared with normal tissue. An absolute log2 fold change (FC) > 1 and an adjusted P value of < 0.05 cutoff was used to defined DELs.

Construction and validation of the eight-lincRNA signature

The DELs of the training set (GSE53625) were exposed to univariate COX regression and 33 DELs related to the OS (P < 0.05) were measured as predictive lincRNAs for LASSO analysis (Table 1). Finally, eight lincRNAs were selected to construct a risk signature for ESCC using multivariate COX regression. The risk score of the signature for OS was identified: risk score = 18.977 × (expression level of AP000487) + 11.606 × (expression level of AC011997) + 0.023 × (expression level of LINC01592) + 3.658 × (expression level of LINC01497) + 24.196 × (expression level of LINC01711) + 3.164 × (expression level of FENDRR) + 0.228 × (expression level of AC087045) + 10.548 × (expression level of AC137770). As above, it was suggested that the lncRNAs in the signature were all risk factors for OS. Their coefficients indicated their impact on OS prediction. For example, the influence of LINC01711 was greatest while that of LINC01592 was least. All patients were separated to high- and low-risk sets based on median risk score in the GEO (Figure 2B) and TCGA (Figure 2C) cohorts. The patients’ status, survival time, and lincRNA expression levels are shown in Figure 2B (GEO) and 2C (TCGA). The survival analysis presented that the OS of low-risk set was better than that of high-risk set in the GEO cohort (hazard ratio, HR = 3.79, 95% confidence interval, 95% CI [2.56-5.62]; P < 0.001) (Figure 3A). The results were consistent in the TCGA cohort (HR = 2.73, 95%CI [1.65-4.53]; P < 0.001) (Figure 3B). The area under the ROC curve (AUC) for 0.5-, 1-, 3-, and 5-year OS were 0.673, 0.734, 0.798, 0.816, 0.795 and 0.777, 0.644, 0.642, 0.649, 0.765 in the TCGA and GEO cohorts, respectively. Together, it was indicated that the signature showed an excellent performance for OS prediction.

Figure 2. Establishment and validation of the eight-lincRNA prognostic signature. (A) The procedure of the establishment of the prognostic signature. (B–C) Correlation between the prognostic signature and the overall survival of patients in the GEO cohort (B) and TCGA (C) cohorts. The distribution of risk scores (upper), survival time (middle) and lincRNA expression levels (below). The black dotted lines represent the median risk score cut-off dividing patients into low- and high-risk groups. The red dots and lines represent the patients in high-risk groups. The green dots and lines represent the patients in low-risk groups.

Figure 3. Kaplan–Meier survival and ROC curves of the eight-lincRNA prognostic signature. (A, B) Kaplan-Meier survival curves of overall survival among risk stratification groups in the GEO (A) and TCGA (B) set. (C, D) ROC curves with calculated AUCs for risk prediction in 0.5-, 1-, 2-, 3-, 5-years in the GEO (C) and TCGA (D) sets.

Subgroup, univariate and multivariate COX regression analysis of the signature

As shown in Figure 4, the risk score identified by prognostic signature assists as a helpful biomarker for OS prediction in different subgroups, including stage I,II (P < 0.0001), stage III (P < 0.0001), age≤60 (P < 0.0001), age > 60 (P < 0.0001), grade1,2 (P < 0.0001), grade 3 (P < 0.0001) in the GEO cohort, and stage I,II (P = 0.012), stage III (P = 0.006), T1,2 (P = 0.013), T3,4 (P = 0.01), age ≤ 60 (P = 0.012), age > 60 (P = 0.004) in the TCGA cohort, respectively.

Figure 4. Confirmation of the signature stratified by different clinical factors in the GEO and TCGA cohorts. Kaplan–Meier survival for OS in subgroups stratified by stage I,II (A) stage III (B) age≤60 (C) age > 60 (D) grade 1,2 (E) grade 3 (F) in the GEO cohort, and stage I,II (G) stage III (H) T1,2(I) T3,4(J) age≤60 (K) age > 60 (L) in the TCGA cohort.

The univariate Cox regression showed that grade, stage, N stage and risk score in the TCGA cohort (grade: P =0.048, stage: P < 0.001; N stage: P < 0.041, risk score: P < 0.001; Figure 5A), and stage, risk score in the GEO cohort (stage: P < 0.001; risk score: P < 0.001; Figure 5C) were predictors for OS. Moreover, multivariate Cox regression analysis confirmed that age (HR = 1.032, 95% CI [1.008–1.056]; P = 0.009; Figure 5B) and risk score (HR = 1.617; 95% CI [1.444–1.811]; P < 0.001; Figure 5B) were significant independent risk factors in the TCGA cohort. Multivariate Cox regression further showed that stage (HR = 2.401, 95% CI [1.606–3.590]; P < 0.001; Figure 5D) and risk score (HR = 1.142; 95% CI [1.063–1.226]; P < 0.001; Figure 5D) was significant independent risk factors in the GEO cohort. These data indicated that the signature was an independent risk factor of ESCC.

Figure 5. Univariate and multivariate Cox regression analyses of clinical factors associated with overall survival. (A-B) Univariate Cox regression analyses of clinical factors associated with overall survival in the TCGA (A) and GEO (C) set. (C-D) Multivariate Cox regression analyses of clinical factors associated with overall survival in the TCGA (B) and GEO (D) sets.

Nomogram construction

Based on the prognostic signature and clinical factors, such as age, grade and stage, a nomogram was constructed (Figure 6A). The calibration curve was used to describe the prediction value of the nomogram and the 45-degree line indicated the actual survival outcomes. The results for predicting 1-, 3- and 5-year OS showed that the nomogram-predicted survival closely matched with the best prediction performance (Figure 6B). The 1-year AUC was 0.734 for nomogram, and 0.611 for age, 0.530 for grade, 0.584 for stage. The 3-year AUC was 0.816 for nomogram, and 0.591 for age, 0.600 for grade, 0.646 for stage. Moreover, the 5-year AUC was 0.795 for nomogram, and 0.589 for age, 0.584 for grade, 0.639 for stage. These findings showed that compared with a single clinical factor, the nomogram combined the signature and clinical factors had great predictive accuracy.

Figure 6. Construction of a nomogram for overall survival prediction in ESCC. (A) The nomogram consists of age, grade, stage and the risk score based on the eight-lncRNA signature. (B) Calibration curves of the nomogram for the estimation of survival rates at 1-, 3-, 5- year. (C) The Kaplan-Meier curves of the risk subgroups stratified by the tertiles of total points derived from the nomogram.

WGCNA

WGCNA was used to develop a gene co-expression system to select biologically meaningful gene modules which are related to the lincRNAs in the signature. We set the cut-off as Person correlation coefficient > 0.6 and P < 0.001 to screen genes co-expression with lincRNAs. Then, differential expressed gene (DEG) analysis was performed among these genes. Figure 7A shows the volcano plot of DEGs. To create a scale-free system, the scale-free topology fit index reached 0.8 by setting the soft threshold power value beta to 5 (Figure 7B). In addition, genes with similar patterns was located in different modules by average linkage clustering (Figure 7C). The minimum cluster size was identified as 30 per module. The dynamic shear method was performed to determine the gene module. We calculated module eigengene (ME) which was the overall gene expression level of corresponding modules and clustered them based on their correlation in order to explore co-expression similarity of all modules (Figure 7D). To analyze the connection of gene modules and clinical characteristics, eigengenes were calculated correlations with clinical factors, such as survival time, status, gender, age, alcohol, grade, stage, T, N, pneumonia. Finally, the robust correlation was identified between the gene significance and the survival time, T stage (Figure 7E). The eight modules were generally separated into two clusters (Figure 7F). Gene significance were also determined to evaluate the correlation between gene expression and survival time (Figure 8A). To further analyze genes in these modules, we found a strong correlation between the gene significance for the survival time and Module Membership in the green module (cor = 0.47, P = 2.8e − 08), the brown module was negatively correlated with the survival time (cor = −0.24, P = 0.0025), the blue module was positively correlated with the survival time (cor = 0.45, P = 2.9e − 11), and the red module was positively correlated with the survival time (cor = 0.3, P = 0.0026) (Figure 8B). The function enrichment analysis was performed to explore the GO database and KEGG pathway in which are involved (Figure 8C). The results indicated that the biological process mainly involved in keratinocyte differentiation, peptide cross-linking, epidermal cell differentiation, and skin development and so on. The results showed that the molecular function was related to peptidase activity, acting on L-amino acid peptides, endopeptidase inhibitor activity, cadherin binding and serine-type peptidase activity. The cell components which was correlated to the results included tertiary granule lumen, intermediate filament cytoskeleton, keratin filament, and intermediate filament. Furthermore, KEGG pathway functional enrichment showed that influenza A, chemical carcinogenesis and drug metabolism were mainly related to the genes in these modules.

Figure 7. WGCNA analysis. (A) Volcano plot of differentially expressed genes associated with lincRNAs in the signature (correlation coefficient < 0.7, P < 0.001). (B) Analysis of the scale-free topology model fit index for various soft-thresholding powers (β) and the mean connectivity for various soft-thresholding powers. In all, 5 was the most fit power value. (C) Dendrogram of the genes and different clinical factors of ESCC (survival time, status, gender, age, alcohol, grade, stage, T, N, Pneumonia). (D) Dendrogram of the gene modules based on a dissimilarity measure. The branches of the cluster dendrogram correspond to the different gene modules. Each piece of the leaves on the cluster dendrogram corresponds to a gene. (E) Module-trait relationships. Heatmap of the correlation between module eigengenes and clinical characteristics of ESCC. (F) Hierarchical clustering and heatmap of the hub gene network.

Figure 8. The correlation between the genes in the modules and survival time. (A) Distribution of mean gene significance and standard deviation with survival time in the modules of ESCC. (B) Scatter plot of module eigengenes in eight modules. (C) GO and KEGG pathway enrichment of eight modules. GO enrichment contains three categories including biological process, cellular component and molecular function.

The same analysis was completed to assess the correlation between gene expression and T stage (Figure 9A). Thus, we found a strong correlation between the gene significance for the T stage and Module Membership in the green module (cor = 0.57, P = 3.3e − 12), the black module was positively correlated with the T stage (cor = − 0.24, P = 2.1e – 8), and the blue module was negatively correlated with the T stage (cor = 0.21, P = 0.003) (Figure 9B). The function enrichment analysis was also performed to explore the GO database and KEGG pathway in which are involved (Figure 9C). The results indicated that the biological process mainly involved in extracellular matrix organization. The results showed that the molecular function was related to metallopeptidase activity, metalloendopeptidase activity and collagen biding. The cell components which was mainly correlated to the results included endoplasmic reticulum lumen. Moreover, KEGG pathway functional enrichment showed that protein digestion and absorption was mainly related to the genes in these modules.

Figure 9. The correlation between the genes in the modules and T stage. (A) Distribution of mean gene significance and standard deviation with survival time in the modules of ESCC. (B) Scatter plot of module eigengenes in eight modules. (C) GO and KEGG pathway enrichment of eight modules. GO enrichment contains three categories including biological process, cellular component and molecular function.

Data collection

The RNA-sequencing and clinical information for ESCC were acquired from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) as the training cohort, and The Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov/) as the validation cohort. GSE53625 from GEO was conducted by GPL18109 (Agilent-038314 CBC Homo sapiens lncRNA + mRNA microarray V2.0) [28]. The clinical features of patients with ESCC in the TCGA and GEO cohorts are presented in Table 2.

Identification of differentially expressed lincRNAs in ESCC

The RNA-sequencing data were normalized with Expectation-Maximization algorithm (log2 transformation). Then 37501 lincRNAs were annotated. The differentially expressed lincRNAs (DELs) were identified between tumor and normal tissues using the R package “Limma” with log2 | fold-change (FC) | > 1 and adjusted P-value < 0.05.

The construction of lncRNAs-based prognostic signature

First, we used univariate and least absolute shrinkage and selection operator (LASSO) COX regression analysis to select the independent risk lncRNAs. Next, multivariate COX regression was used to identify corresponding coefficients of ESCC prognostic signature using R package “glment”, “survminer” and “survival”. The risk score of every patient from the TCGA and GEO cohorts were calculated based on the signature. All samples were randomly separated to high- and low-risk sets with the median score as cut-off value. Survival analysis for each set was evaluated by the Kaplan–Meier curve and log-rank test. The receiver operating characteristic (ROC) curve and the area under the curve (AUC) were drawn using R package “survivalROC”.

Weighted gene co-expression network analysis (WGCNA)

WGCNA is a useful tool to establish the co-expression network between gene pattern and clinical factors using R package “WGCNA” (Version: 1.68). The procedure of WGCNA included identifying gene expression similarity matrix, adjacency matrix and co-expression network. ScaleFree plot was used for evaluating whether the network exhibits a scale free topology. The power value of soft threshold of adjacency matrix was identified as 5 to meet scale-free topology criterion. The hierarchically clustering analysis based on average-linkage were originated from the Dynamic using Tree Cut method for Branch Cutting (deep-split = 2, cut height = 0.4, minimum cluster size = 30). The association between modules and variables was identified to select relevant module.

The nomogram establishing

A concise nomogram of predicting the survival of ESCC was established using the R package “rms”, “Hmisc”, “lattice”, “Formula”, and “foreign”. The concordance index (C-index) was performed to assess prediction capability. The patients with ESCC were separated to diverse risk clusters along with their scores. All analyses were completed with R software. The P-value < 0.05 was considered statistically significant.

Data availability

All datasets generated for this study are included in the manuscript.

Ethics statement

As the data (TCGA and GEO datasets) are publicly available, no ethical approval is required.