Exploring the Potential of Exosome-Related LncRNA Pairs as Predictors for Immune Microenvironment, Survival Outcome, and Microbiotain Landscape in Esophageal Squamous Cell Carcinoma

The RNA sequencing data and clinical information of ESCC patients were obtained from the TCGA database (80 tumor and 11 normal samples) and the GEO database (GSE53624↗, 119 tumor samples). The number of transcriptional samples (Stage I-II: 53 samples; Stage III-IV 64samples) is balanced in the GEO dataset. In each dataset, duplicate sequencing samples from the same patient were excluded, and patients lacking complete follow-up information and with zero survival days were excluded. For TCGA-ESCC cohort, convert RNA-seq data in FPKM format to TPM. The GEO cohort is annotated by the Affymetrix platform. The ComBat algorithm in the “sva” package is used for batch effect elimination (meta queue) for the TCGA and GEO databases. The GENCODE database (GRCh38 version) is used for lncRNA and mRNA annotation, taking the intersection of the two platforms. 121 exosome-related genes (ER-mRNAs) were downloaded from the ExoBCD database (https://exobcd.liumwei.org/↗), and finally 94 ER-mRNAs could be annotated in the Meta cohort. The above ER-mRNA was used to screen exosome-related lncRNAs (ER-lncRNAs) by co-expression strategy with correlation coefficient > 0.4 and p-value < 0.001 as the screening threshold. The “limma” package was used for differential expression analysis between ER-lncRNAs and ER-mRNAs with thresholds set to logFC >1 and FDR <0.05 Real.

Cyclic pairing was performed on the differential ER-lncRNAs, assuming C=lncRNA-A + lncRNA-B. If the expression level of lncRNA-A is higher than that of lncRNA-B, the result of C is 1, otherwise the result of C is 0, and finally a new matrix (0 or 1) containing C is generated. If the expression level of a certain lncRNA in ER-lncRNA Pairs is 0 or 1, the relationship between Pairs and prognosis is not considered. A valid match was considered when the number of ER-lncRNA Pairs expressed as 0 or 1 accounted for more than 20% of the total number of pairs. The flow chart is shown in Figure 1A.

According to the expression profiles of ER-lncRNA, consensus clustering was performed using the “ConsensusClusterPlus” package in R by setting the number of groups to 9, 100 replicates, and the pltem=0.8. The optimal cluster number was calculated using the consensus matrix and cumulative distribution function (CDF). The differences in the overall survival (OS) between different clusters were estimated using the Kaplan-Meier method. Comparisons of the distribution of categorical data among the clusters were done using the chi-squared test.

The GSVA method was used to assess differences in biological pathways between subtypes. Differential pathways between different subtypes were analyzed using the “limma” package, with c2.cp.kegg.v7.0.symbols.gmt as the reference gene set and FDR < 0.05 as the screening threshold.

Univariate Cox regression analysis was performed for all ERlncRNA pairs (P < 0.001), followed by least absolute shrinkage and selection operator (LASSO) Cox regression analysis (10-fold cross-validation). LASSO regression was run for 1000 cycles, and for each cycle, 1000 random times were set to remove redundant ER-lncRNA pairs. We included the above ER-lncRNA pairs in the Cox regression risk model. When the AIC value was minimal, the calculation was stopped and the model was also used as the best model. The risk score was calculated as follows:

The cut-off values were determined by evaluating the 1-year ROC curves in the training set (GEO cohort) and the validation set (TCGA cohort). Patients were divided into high-risk and low-risk groups according to the cut-off values. Subsequently, Kaplan-Meier curves and ROC curves were applied to assess the prognostic role of the model. To verify the clinical application value of the constructed model, we analyzed the association between the model-based risk score and clinicopathological features based on the GEO database. Additionally, univariate and multifactor Cox regression analyses were performed to determine whether the model could be used as an independent predictor of clinical prognosis.

The ssGSEA algorithm was used to calculate the scores of 29 immune cells and immune function in ESCC samples, and the “pheatmap” package was used for heatmap visualization. The ESTIMATE algorithm identifies specific features associated with stromal cell and immune cell infiltration. Differences between different molecular subtypes and risk groups were compared using the Kruskal Wallis test. We used TIMER, XCELL, MCPcounter, EPIC algorithms to analyze the relationship between risk scores and immune cell characteristics. Additionally, we calculated enrichment scores for immune microenvironment-related pathways in each ESCC sample using the “IOBR” package.

IC50 was calculated using the “prophetic” package in R software, and chemotherapeutic drugs were obtained from the Genome of Drug Sensitivity in Cancer (GDSC) database.

A total of 10 tumor tissue samples and nearby normal esophageal tissue samples were obtained from ESCC patients who underwent tumor resection. All tissue samples were collected from the Thoracic Surgery Department of the Second Hospital of Hebei Medical University with the approval by the Medical Ethics Committee of the hospital. The specific experimental protocol comes from our previous research (21, 22).

All the statistical analyses were made with the R software (v.4.0.1). The above section has described detailed statistical approaches for transcriptome data processing. A p-value less than 0.05 was of statistical significance.

We identified 41 differentially expressed ER-mRNAs in normal and tumor samples in the meta cohort following the details described in the methods (Figures 1B, C). Subsequently, co-expression analysis was performed between differentially expressed ER-mRNAs and lncRNAs in tumor samples from the meta cohort, and a total of 620 ER-lncRNAs were identified. Ultimately, secondary differential analysis was performed in different samples, and 142 were defined as differentially expressed ER-lncRNAs (Figure 1D), of which 109 were upregulated and 33 were downregulated (Figure 1E).

According to the method described, we identified 3459 valid ER-lncRNA pairs in the meta cohort (Supplementary File 1). These ER-lncRNA pairs were further used for clustering analysis. The cluster effect was best when ESCC patients were clustered into three subtypes, and the subtype internal consistency and stability were good (Figures 2A–C). PCA analysis revealed significant heterogeneity in molecular typing based on ER-lncRNA pairs (Figure 2D). Survival analysis showed that subtype B had the worst prognosis of these three molecular subtypes, whereas subtype A had the best prognosis (Figure 2E). The ESTIMATE algorithm was used to assess stromal score, immune score and tumor purity for the three subtypes. As shown in Figures 2F–H, subtype C had the highest stromal and immune score and the lowest tumor purity score, while subtype A had the lowest stromal and immune score and the highest tumor purity score.

Of the three molecular subtypes, subtype C showed the highest enrichment and activity in several immune cells, pathways or functions (Figure 3A). Subsequently, we examined the expression of four immune checkpoint genes (i.e., PDCD1, CTLA4, HAVCR2, and LAG3) and found that the mRNA levels of most of the immune checkpoints were higher in subtype C than in other subtypes, which may suggest that this subtype may benefit more from immunotherapy (Figure 3B). In addition, we used the pRRophetic algorithm to evaluate the therapeutic effect of chemotherapy drugs. Among some classical chemotherapy drugs, subtype B was more sensitive to cisplatin, doxorubicin, gemcitabine, paclitaxel (Figure 3C), suggesting that molecular subtype-based grouping may play a guiding role in chemotherapy.

To investigate the causes of the different survival states and immune landscapes, we used GSVA to study the biological processes between the different subtypes. The results showed that subtype A had more activated pathways compared with subtype B, such as KEGG_NOD_LIKE_RECEPTOR_SIGNALING_PATHWAY, KEGG_SPHINGOLIPID_METABOLISM, KEGG_ETHER_LIPID_METABOLISM, etc. (Figure 4A). Compared with subtype C, subtype A has more activated pathways, such as KEGG_TIGHT_JUNCTION, KEGG_TASTE_TRANSDUCTION, etc. (Figure 4B). Compared with subtype B, subtype C has more activated pathways, such as KEGG_GLYCOSAMINOGLYCAN_BIOSYNTHESIS_CHONDROITIN_SULFATE,KEGG_EPITHELIAL_CELL_SIGNALING_IN_HELICOBACTER_PYLORI_INFECTIO, etc. (Supplementary Figure 1)

We performed univariate Cox regression analysis in the GEO cohort (P < 0.001) and identified 8 ER-lncRNA pairs with strong prognostic significance (Figure 5A). Eight ER-lncRNA pairs (AP000487.1↗ |AC011483.1↗, AC127526.2↗ |C8orf49, ADAMTS9-AS2|PLA2G4E-AS1, AC104653.1↗ |PLA2G4E-AS1, AL356056.2↗ |IFNG-AS1, PLA2G4E-AS1|LINC00460, AC116903.2↗ |AC082651.3↗ 和LINC00958|AP002336.2↗) were included in the Cox regression model through LASSO regression analysis. Next, we calculated the AUC values of the risk-score model based on 8 pairs in the GEO cohort, where the cut-off values for 1 (Figure 5B), 2 (Figure 5C), 3 (Figure 5D), and 4 years (Figure 5E) were 1.135, 1.135, 0.504, and 0.504, respectively. Additionally, based on the same risk-score model, the validation set demonstrated superior survival prediction ability (Figures 5F, G). In addition, we grouped patients by their 1-year cut-off values in each cohort. Kaplan-Meier analysis showed that in the GEO cohort (Figure 5H) and the TCGA cohort (Figure 5I), patients in the low-risk group survived longer than those in the high-risk group (p < 0.001).

The risk scores and survival rates for each case in the GEO cohort are shown in Figures 6A, B. These data suggest that patients in the low-risk group had better clinical outcomes than those in the high-risk group. In addition, ROC curve analysis showed that for median survival time, ER-lncRNA pairs–based risk stratification had better predictive indicative value than traditional clinical indicators (Figure 6C). The results of the DCA analysis also showed that ER-lncRNA pairs–based risk stratification had a better level of benefit relative to traditional clinical indicators (Figure 6D).

To estimate the correlation between risk-score model and clinicopathological characteristics of ESCC patients, we performed chi-square test in the GEO cohort, which showed that grade, N stage, and stage were significantly related to the risk-score model (Figure 7A). Furthermore, according to the Wilcoxon rank-sum test, the risk scores of ESCC patients were significantly related to clinical stage, status of differentiation, and status of lymph node metastasis (Figures 7B–D). Additionally, in the GEO cohort, we demonstrated that the risk score can be used as an independent prognostic factor for ESCC patients (Figures 7E, F). Consistent results were also exhibited in the TCGA cohort (Supplementary Figure 2).

To further explore the relationship between risk assessment model and the tumor immune microenvironment, we analyzed the composition of immune cells in different risk subgroups in the meta cohort by applying multiple deconvolution algorithms. The results showed that the high-risk group had a higher landscape of immune cell infiltration (Figure 8A), as well as a more active TME (Supplementary Figure 3), compared to the low-risk group. In addition, in immune checkpoint mRNA expression, we found that high-risk patients had higher expression of PDCD1, CTLA4, and HAVCR2 (Figure 8B). Meanwhile, high-risk patients also had smaller IC50 values among most chemotherapeutic agents (Figure 8C), suggesting that the model could serve as a potential predictor of chemosensitivity.

The Cancer Microbiome Atlas (TCMA) (23), which reveals the distribution of microbiota in TCGA sequenced samples, enables a systematically matched microbial-host multi-omics analysis. Our secondary analysis of the ESCC data among them founded that the ESCC environment contains the following five main microbial phyla: Firmicutes, Bacteroidetes, Proteobacteria, Fusobacteria, and Actinobacteria (Figures 9A, B). Interestingly, Actinobacteria was significantly upregulated in high-risk patients (Figure 9C), which may suggest us a large role of Actinobacteria in the development of ESCC.

We split the 8-pairs in the risk-score model and finally got 14 ER-lncRNAs. According to the method in the literature (24), with survival and death as binary dependent variables, based on the expression of 14 ER-lncRNAs, a variety of machine learning algorithms (DT, GBM, LF, NNET, RF and XGBOOST) were used for prediction. It was found that the AUCs of XGBOOST and RF were the largest (Figure 10A). In the RF model, the five most important ER-lncRNAs were: AC082651.3↗, AP000487.1↗, PLA2G4E-AS1, C8orf49 and AL356056.2↗ (Figure 10B). Subsequently, the expression levels of the above lncRNAs were analyzed by combining the GTEx and TCGA databases, in which only C8orf49 was upregulated in normal samples, and the rest of lncRNAs were upregulated in tumor samples (Figure 10C). In addition, qRT-PCR analysis based on clinical samples from our hospital showed a high degree of consistency (Figure 10D).

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/. 1