Identification of m6A/m5C-related lncRNA signature for prediction of prognosis and immunotherapy efficacy in esophageal squamous cell carcinoma – Scientific Reports

Screening m6A/m5C-related lncRNAs

Using the defined criteria (|r|> 0.3 and P < 0.05), we identified a total of 2091 m5C gene-related lncRNAs and 2366 m6A gene-related lncRNAs within ESCC samples compared with normal samples. The top 30 lncRNAs related to m5C and m6A genes were listed in the correlation matrix (Fig. 1A,B). Based on the sequencing data in TCGA and GEO dataset, the lncRNA included in both two datasets were selected. Next, the overlapped lncRNAs between m5C and m6A gene-related lncRNAs were obtained as candidate molecules. Finally, we screened a total of 606 m6A/m5C-related lncRNAs for further analysis (Fig. 1C).

Figure 1
figure 1

Identification of m6A/m5C-related lncRNAs in ESCC. (A,B) Top 30 m5C-related lncRNAs and top 30 m6A-related lncRNAs. (C) Overlapping lncRNAs between m5C- and m6A-related lncRNAs. (D) Functional enrichment analysis revealing KEGG pathways enriched by the lncRNA target genes. (EG) Enrichment of GO biological process, cellular component, and molecular function. (H) Cellular localization of m6A/m5C-lncRNAs.

For these lncRNAs, we predicted their potential target genes by examining coding genes located 10 km upstream and downstream of each lncRNA, obtaining 758 target genes. We performed functional analysis, and found these m6A/m5C-lncRNA target genes were significantly enriched in the olfactory transduction pathway (Fig. 1D). Moreover, these genes demonstrated relevance to multiple biological process, cellular component, and molecular function, such as metabolism-related process, fibrinogen complex, chromatoid body, olfactory receptor activity and odorant binding (Fig. 1E–G).

In addition, we utilized the LncSLdb/lncATLAS database to predict the subcellular localization of lncRNAs. Our findings indicated that 61.55% of m6A/m5C-lncRNAs were localized in the ribosome, while 33.17% of lncRNAs were situated in the cytoplasm. Moreover, 4.79% of m6A/m5C-lncRNAs were found in exosomes, and 0.5% were identified in the nucleus (Fig. 1H).

Molecular classification of ESCC case based on m6A/m5C-lncRNAs expression

Based on the expression differences of m6A/m5C-lncRNAs, we performed unsupervised classification of the ESCC cases in the TCGA dataset, ultimately identifying three distinct clusters (Fig. 2A,B). Survival analysis demonstrated significant variations in survival outcomes among these three clusters (Fig. 2C). Specifically, patients in cluster 1 exhibited superior survival outcomes compared to those in the other two clusters.

Figure 2
figure 2

Identifying three clusters of ESCC cases by consensus clustering in the TCGA set. (A) Consensus clustering analysis and relative change in the area under the CDF curve. (B) Heat map of consensus matrix. (C) Prognosis analysis among the three clusters. (D) Heatmap of the clinical features among the three clusters in ESCC cases. (E) GSVA analysis displaying the key pathways enriched by 50 hallmark gene sets across three clusters. (FG) Box plots illustrating differences in immunity landscape and immune checkpoint gene expression among three clusters.

We further investigated the differences in various clinical features among the subtypes within the TCGA-ESCC dataset. The matrix results indicated that the m6A/m5C-lncRNAs expression in cluster 3 was relatively lower. In contrast, cluster 2 displayed moderate expression, while cluster 1 exhibited the highest expression levels of m6A/m5C-lncRNAs. Additionally, differences were observed in clinical pathological parameters, including age, grade, immune subtype, and gastric classification, across the three clusters (Fig. 2D).

Notably, significant differences were identified among the three clusters in the hallmark gene set. For instance, samples in cluster 1 displayed relatively lower GSVA enrichment scores across the 50 hallmark gene sets, while cluster 3 exhibited high scores. A total of 32 hallmark gene sets exhibited significant differences between the three clusters (Fig. 2E).

We analyzed immune landscape differences among the three distinct clusters using the xCell, CIBERSORT, and ssGSEA methods for the TCGA-ESCC samples. Our findings indicated that the infiltration ratio of CD4 + memory T cells, CD8 + T cells, and other T cell-related cell types was significantly lower in cluster 3, as evident in the xCELL results (Fig. 2F,G). Furthermore, a total of 17 immune checkpoint genes exhibited significant expression difference among the three clusters. Notably, most of these genes were markedly upregulated in cluster 1 compared to cluster 3 (Fig. 2G).

Construction of m6A/m5C-lncRNA-based prognostic signature

To identify the prognostic lncRNAs, we conducted a cox regression analysis based on the TCGA-ESCC samples. There were 17 m6A/m5C-lncRNAs with significant prognostic differences (Table 2). We classified the ESCC cases into higher and lower expression groups using the median gene expression value as the criterion. Survival analysis results demonstrated significant prognostic differences among 8 of these lncRNAs between the high- and low-expression groups (Fig. 3).

Table 2 Univariate cox regression analysis showing the m6A/m5C-related prognostic lncRNAs.
Figure 3
figure 3

Survival analysis results of prognostic lncRNAs between high- and low- expression groups of ESCC samples.

According to lasso-cox analysis, we finally obtained 10 prognostic lncRNAs, including LINC00847, LINC00942, TTTY15, LINC01602, LINC01310, LINC00898, EGFR-AS1, MIF-AS1, LINC00993, and LINC01415. A ten lncRNAs-based RiskScore model was generated for survival prediction (Fig. 4A–C).

Figure 4
figure 4

Construction of lncRNAs signature-based RiskScroe system for prognosis prediction. (A,B) Construction of the RiskScore model using the lasso method. (C) Lasso coefficients of the 10 prognostic genes. (D) Survival analysis of ESCC patients in high- and low- RiskScore cohorts. (E) ROC curve of the RiskScore at 1-, 2- and 3-year follow-up. (FH) Distribution of RiskScore, status, and expression of m6A/m5C-lncRNAs in the TCGA.

In the training set, we observed that patients in the high-RiskScore subtype had a worse prognosis than those in the low-RiskScore group (Fig. 4D,F–H). Furthermore, ROC curve analysis indicated that this model exhibited promising predictive ability for survival in TCGA-ESCC samples (AUC value at 1-, 3-, and 5- years were respectively 0.805, 0.921 and 0.903; Fig. 4E).

Validation of the RiskScore model

The predictive ability was also evaluated in the validation set GSE53622. There was a significant prognostic difference between the two subgroups (Fig. 5A,C–E). ROC curve confirmed the superior predictive power of this RiskScore model in the validation sets (AUC values at 1-, 2-, and 3- years were respectively 0.545, 0.583, 0.590; Fig. 5B).

Figure 5
figure 5

Prognostic value of RiskScore models in an external validation set. (A) Survival analysis of two subgroups via the GSE53622 dataset. (B) ROC curve analysis of the RiskScore model in the validation set. (CE) Distribution of RiskScore, status, and the expression level of lncRNAs in the GSE53622 dataset.

Next, the prognostic efficiency of the RiskScore system in various clinical feature subgroups was confirmed on the TCGA-ESCC samples. The results demonstrated significant survival differences between the high- and low- RiskScore groups in different clinical characteristics, age, grade, stage, immune subtype, and SCNA level (Fig. 6A).

Figure 6
figure 6

Prognostic efficiency prediction of m6A/m5C-lncRNA signatures. (A) Predictive value of RiskScore in clinicopathological subgroups. (B) RiskScore assessments across different clinicopathological groups, including survival status, age, grade, gender, stage, somatic copy number alterations (SCNA) levels, and gastric classification. (C) Sankey diagram visualized the correlation ship of clusters, RiskScore group and survival status. (D) Risk score assessments between different cluster cohorts.

Furthermore, ESCC patients in the TCGA dataset could be classified into different subgroups based on clinical features, such as stage, grade, and age. RiskScore assessments were performed between different clinical feature groups and the previously identified three clusters. Notably, patients in the alive and female subgroups presented a low RiskScore, indicating a significant correlation between RiskScore and survival status and gender. Moreover, a significant difference in RiskScore was observed between the G1 and G2 grade subgroup (P < 0.05, Fig. 6B). Patients in cluster 1 displayed the lowest RiskScore than the other two clusters (P < 0.05, Fig. 6D). Sankey diagram showed the consistency between cluster and RiskScore for survival status prediction, which exhibited a consistent result to the RiskScore assessments (Fig. 6C,D).

Identification of RiskScore as an independent prognostic factor

Utilizing survival information from the TCGA-ESCC dataset, we conducted a cox regression analysis, and the results revealed that RiskScore emerged as an independent prognostic factor for survival prediction (Fig. 7A,B; P < 0.05, P < 0.01).

Figure 7
figure 7

Cox regression analysis and nomogram model construction for m6A/m5C-lncRNAs prognostic signature. (A,B) Cox regression analysis for the survival time in TCGA cohorts. (C) The nomogram model predicting survival outcomes of ESCC. (DG) ROC and calibration curves for predicting the survival time in the TCGA cohort.

We constructed a nomogram model that integrated clinical characteristics and RiskScore for predicting survival in the TCGA (Fig. 7C). To gauge the predictive accuracy, we employed ROC curves, revealing AUC values of 0.518 and 0.648 for age and stage, respectively (Fig. 7G). The nomogram prediction and RiskScore possessed higher AUC values (0.886, 0.818) than other clinically relevant prognostic factors in the TCGA cohort. It is worth noting that due to limited data, we refrain from drawing generalized conclusions beyond 3 years. The calibration curves demonstrated consistency with the standard curve in predictive outcome at 0.5, 1, and 2 years (Fig. 7D–F). These findings revealed nomogram model exhibited reliable prognostic predictive power compared with other clinical characteristics in TCGA cohorts.

Association of RiskScore and several highly trustworthy indices

Cancer progression encompasses the transition towards a dedifferentiated oncogenic phenotype and the acquisition of stem cell-like features. A previous study employed the innovative logistic regression OCLR algorithm to introduce two stemness indices for quantifying tumor microenvironment stemness22: mRNA expression and the methylated DNA-based stemness index (mRNAsi and mDNAsi). These indices reflect the gene expression and epigenetic characteristics of stem cells. Epigenetic regulation-based mRNAsi and mDNAsi (EREG-mRNAsi and EREG-mDNAsi) were generated via reconstructing the gene regulatory network from methylation and transcriptome data22. Additionally, single nucleotide variants (SNVs) represent the most prevalent genomic mutations in tumor cells, leading to the production of variant peptides distinct from wild-type proteins, which are subsequently presented by MHC-I23.

In this context, we explored the correlation between RiskScore and several highly reputable indices within the TCGA-ESCC cohort. These indices encompassed stemness indices (mRNAsi, mDNAsi, EREG-mDNAsi, DMPsi), tumor mutational burden (TMB), stromal score, tumor purity, immune score, and SNV-neoantigen. Our findings unveiled a negative correlation between the m6A/m5C-lncRNA signature and mDNAsi, EREG-mDNAsi, DMPsi, ENHsi, TMB, and SNV-neoantigen. Further analysis revealed significant statistical differences between the two groups in mDNAsi, DMPsi, and EREG-mDNAsi. Particularly noteworthy was that ESCCs in the high-RiskScore groups exhibited lower mRNAsi, DMPsi, and EREG-mDNAsi compared to the low-RiskScore group (Fig. 8, P < 0.05).

Figure 8
figure 8

Correlation analysis between RiskScore and several clinical trustworthy indices. *P < 0.05; **P < 0.01; ***P < 0.001.

Correlation between lncRNA signature and immune cell proportion

We evaluated the immune cell ratio in the TCGA-ESCC dataset between the two risk groups. Based on xCell results, we observed a significantly increased ratio of CD4 + T cells in the low-RiskScore group compared to the high-RiskScore group (Fig. 9A, P < 0.05). Conversely, the cibersort results indicated a decreased proportion of activated NK cells and macrophage M2 in the low-RiskScore group (Fig. 9A, P < 0.05). However, ssGSEA analysis revealed no statistical difference in infiltrated immune cell levels between the high and low-RiskScore groups (Fig. 9A, P < 0.05).

Figure 9
figure 9

Association of RiskScore, infiltrated immune cell and immunotherapy response via TCGA-ESCC datasets. (A,B) Differences in immune cell infiltration and immune checkpoint gene expression between two subgroups. *P < 0.05; **P < 0.01; ***P < 0.001. (C) Differences in immunotherapy efficacy between two risk subgroups in TCGA-ESCC cohorts. (D) Immunotherapy response between high and low-RiskScore groups in the IMvigor210 datasets.

Additionally, we analyzed the expression of immune checkpoint genes between the high and low RiskScore groups. Our findings revealed that five immune checkpoint genes exhibited significantly abnormal expression between the two groups, with their expression levels in the low-risk group being notably higher than those in the high-RiskScore group (Fig. 9B).

Moreover, we assessed the predictive value of RiskScore in immunotherapy efficiency using the TIDE method on TCGA-ESCC data. The results indicated that the low-RiskScore group was more likely to respond to immunotherapy compared with high- RiskScore patients (P = 0.025, Fig. 9C). Similarly, we analyzed the proportions of patients who responded to immunotherapy in the high and low-RiskScore groups using the IMvigor210 dataset. Our findings revealed a significant difference in the percentage of drug response patients between the two subgroups (P = 0.032, Fig. 9D), underscoring the reliable predictive performance of the RiskScore model in immunotherapy response.