Development of biomarker signatures associated with anoikis to predict prognosis in patients with esophageal cancer: An observational study

Gene expression profiles, somatic mutation data, copy-number variation data, and clinical details (such as sex, age, survival duration, survival status, and tumor stage) were sourced from The Cancer Genome Atlas database. A total of 183 esophageal cancer samples and 13 samples from healthy subjects were downloaded. Thedataset, comprising 119 tumor samples with clinical data, was sourced from the gene expression omnibus repository.Using the R package “limma,” we identified the differentially expressed genes (DEGs) in TCGA-ESCA. The adjusted-values were evaluated with the false discovery rate correction method, and DEGs were determined by setting the cutoff at |log2FC|>1 and an adjusted-value of ≤ 0.05. ARGs were extracted from GeneCards (. genecards.org). A total of 652 ARGs were identified and listed (Supplemental Digital Content,). The “VennDiagram” package in R was employed to display the intersecting genes between DEGs and ARGs in ESCA. Afterward, we examined the relationship of ARGs with somatic mutation frequency, genetic location, and copy-number variation of the ARGs. Since the research did not include experiments on cells, tissues, or animals, there were no ethical considerations necessary. GSE53624 https://www↗ http://links.lww.com/MD/N605↗ [] ²⁴ p P

Patients with ESCA were grouped into clusters using principal component analysis (PCA). The R software package, “ConsensusClusterPlus”was used to identify the unsupervised clusters of TCGA-ESCA datasets. To validate the clustering, the techniques PCA, T-distributed stochastic neighbor embedding, and uniform manifold approximation and projection were employed using the R packages “broom,” “Rtsne,” and “UMAP,” respectively. The R packages “Survival” and “‘survminer”’ were employed to examine Kaplan–Meier survival curves for various clusters, while “ComplexHeatmap” was utilized to display the gene expression patterns of these clusters. [] ²⁵

Differentially expressed ARGs were collected based on their interaction with the DEGs. Using the R software package “SVA,” the 2 data sets (ESCA and) are merged to avoid batch effects. To further identify prognostic genes associated with ARG, a univariate Cox regression analysis was conducted. Later, least absolute shrinkage and selection operator and multivariate regression analyses were employed to remove genes causing overfitting. After calculating the risk scores according to ARGs, all cases were classified according to their median risk scores. The method for calculating the risk factor involved the following formula: risk factor = (ARGexpression × Coe) + (ARGexpression × Coe) +…+(ARGexpression × Coe), where Coerepresents the correlation coefficient of prognostic genes and ARGrepresents the expression of prognostic genes. Patients diagnosed with ESCA were categorized into groups based on high- and low-risk. The analysis utilized the R packages “survival,” “forest plot” and “glmnet,” maintaining a significance threshold of.05. GSE53624 1 1 2 2 n n n n f f f f P ≤

Examination of the ESCA cohort indicated genetic variations between the high-risk and low-risk groups. A GO analysis was performed to identify the potential pathways related to these differential genes. The GO analysis encompassed 3 categories: molecular functions, biological processes, and cell composition. The “ClusteProfiler” R packagewas utilized to perform GSEA in order to investigate the distinctions between high-risk and low-risk groups. A-value of 0.05 or less was considered statistically significant. [] ²⁶ P

Kaplan–Meier analysis was employed to assess the OS across various clusters of high-risk and low-risk categories. Receiver operating characteristic analysis was employed to assess the accuracy and reliability of prognostic predictions. To assess if ARG characteristics were standalone indicators of ESCA prognosis, both univariate and multivariate regression analyses were conducted on the clinicopathological attributes and ARG signatures. A cumulative risk curve was utilized to compare the outcomes of high-risk and low-risk groups. A clinicopathological nomogram based on ARG was developed to estimate overall survival at 1, 3, and 5 years. Calibration graphs were employed to evaluate the predictive accuracy of the model and to ascertain its capability in nomogram identification. [] ²⁷

ESTIMATEScorewas utilized to examine the variation in immune cell infiltration between groups at low and high risk. ESTIMATEScore is useful for assessing ImmuneScores, which indicate the abundance of immune cells, and stromal scores, which reflect the quantity of stromal cells. In all ESCA mRNA samples, expression ratio matrixes were obtained using the cell-type identification by estimating relative subsets of RNA transcripts function in the R software. Based on the ARGs median expression levels, ESCA samples were categorized into 2 groups, and the importance of immune cell variations between the high and low expression groups was assessed. A-value of 0.05 or less was considered statistically significant. [] ²⁸ P

In order to explore the connection between high- and low-risk groups and their sensitivity to chemotherapy drugs, data on drug responses and targeting pathways were gathered from the genomics of drug sensitivity in cancer. The Spearman correlation analysis of drug information and risk score, the relationship between risk score and drug sensitivity, and the association with the signaling pathway targeted by these drugs were analyzed. The R package “pRRophetic” was utilized to determine the semi-maximum inhibitory concentration (IC50) based on ARG traits, demonstrating their predictive capacity.The IC50 value can be used to measure the apoptosis-inducing ability of a drug. Generally, the stronger the induction ability, the lower the IC50 value. The IC50 value provides information on the resistance of a cell to a drug. [] ²⁹

The R software (4.2.2) was utilized to conduct these analyses. Predictive elements were assessed through Cox regression analysis and the K-M analysis. We carried out survival analysis through K-M estimation and compared groups using the Wilcoxon rank-sum test. A-value of 0.05 or less was considered statistically significant. P

A total of 160 differentially expressed ARGs were obtained from the TCGA-ESCA database, including 18 downregulated and 142 upregulated genes (Figs.A,B). There were 18,687 genes after the 2 datasets (TCGA-ESCA and) with esophageal cancer were merged without the batch effect. Univariate Cox analyses were performed to identify 160 ARGs in the merged datasets. 1 1 GSE53624

To develop an ARG signature, 11 prognostic ARGs were pinpointed through univariate analysis (Fig.C). The network relationships of the ARGs are shown in FigureD. Additionally, the copy-number variation frequencies of the 11 genes were calculated (Fig.E andF). Subsequently, the prognostic ARGs were overmatched using least absolute shrinkage and selection operator regression (Figs.A,B). Finally, 10 ARGs were screened. A multivariate Cox regression analysis was conducted to identify the 5-gene signatures. The calculation process is as follows: risk factor = (MYBBP1A × (−0.4940)) + (PLAU × 0.3281) + (BUB3 × 0.4077) + (HOTAIR × 0.3099) + (EHMT2 × (−0.3481)) (Figs.A,B). 1 1 1 1 2 2 2 2

This equation was determined to calculate the risk score for every ESCA patient. Based on their scores, patients were categorized into groups of high- and low-risk (Fig.C). A significant difference in survival was observed between the 2 groups (Fig.D). The analysis of the ARG signature revealed that the AUC values of the receiver operating characteristic curves for 1-, 3-, and 5-year OS were 0.639, 0.685, and 0.706, respectively (Fig.E). 2 2 2

The PCA analysis was performed to obtain different clusters of ARGs in patients with ESCA. The results showed that k = 3 reached the clustering stability standard (Figs.A–C). As a result, individuals with ESCA can be categorized into 3 clusters: clusters A, B, and C (Fig.D). Kaplan–Meier survival analysis of the 3 clusters found significant prognostic differences in ESCA patients (< .05) (Fig.E). Over time, the survival rate of the patients showed a downward trend, with the prognosis of cluster C being better than that of clusters A and B. FigureF illustrates the distinct clinicopathological features and overall survival rates in ESCA clusters A, B, and C. The results of the ARG cluster difference analysis are shown in FigureG. The results indicated differences between clusters A and B, and clusters B and C (Figs.H,I). 3 3 3 3 3 3 3 3 P

The “ClusterProfiler” package was utilized for GO enrichment analysis to pinpoint genes with differential expression between high-risk and low-risk groups (Fig.A). BPs included cellular responses to starvation, regulation of histone modifications, and responses to starvation. Cell composition included the nucleocytoplasmic transport, peptidase inhibitor, and serine-type endopeptidase complexes. Molecular functions included histone methyltransferase activity to H3-K9, C2H2 zinc finger domain binding, and histone-lysine N-methyltransferase activity. 4

The BPs involved in esophageal cancer with different levels of ARG expression were investigated using GSEA. When comparing the high-risk and low-risk groups, the high-risk group predominantly exhibited signaling pathways related to cytokine-cytokine receptor interaction, graft-versus-host disease, and hematopoietic cell lineage (Fig.B). The signaling pathways linked to low-risk included arachidonic acid metabolism, drug metabolism, cytochrome P450, and fatty acid metabolism (Fig.C). 4 4

To delve deeper into the predictive impact of ARG signatures, an independent analysis was conducted using ARG characteristics and clinical variables from the TCGA dataset. The findings indicated that variables like ARG signature scores and lymphatic node involvement were significant (< .05).This evaluation revealed that the ARG signature risk score continued to be a significant prognostic factor, implying that the suggested anoikis-related signature (ARS) might operate independently of other clinical variables (Fig.A). To create a nomogram model tailored to clinical practice requirements, patient data including sex, age, TN stage, and ARG signatures for those with ESCA were collectively analyzed (Fig.B). Calibration curves (Fig.C) and cumulative hazards (Fig.B). Calibration curves (Fig.C) and cumulative hazards (Fig.D) show a strong correlation between the predicted and OS rates of patients with ESCA. FiguresE–G illustrate the decision curve analysis conducted for the overall survival rates at 1, 3, and 5 years. P 5 5 5 5 5 5 5 5

A clear distinction exists between patients at low and high risk regarding TME, a crucial marker of tumor biological behavior in the ESCA group. Based on the ESTIMATE analysis, low-risk patients had lower ImuneScores, StromalScores, and ESTIMATEScores compared to that of high-risk patients (Fig.A). 6

The cell-type identification by estimating relative subsets of RNA transcripts algorithm was employed to examine the relationship between ARS and immune cell infiltration by analyzing the proportions of 22 different immune cell types (Figs.B–E). Figureb presents a correlation heatmap illustrating the immunotherapy response results for both high-risk and low-risk groups. The results of the immune cell correlation analysis are presented in Fig.C. The analysis of immune correlations identified a link between genes related to immunity and the risk scores of patients (Fig.D). FigureE illustrates the variations in immune response between the high-risk and low-risk groups. A notable statistical disparity was observed between M2 macrophages and inactive mast cells, with their scores showing a negative correlation with patients’ risk factors (Figs.F andG). 6 6 6 6 6 6 6 6

Ultimately, the potential for drug resistance was examined between ESCA patients at high and low risk. In patients with ESCA, IC50 values were correlated with the IC50 scores of targeted and chemotherapy drugs (Fig.). The IC50 values of 16 drugs (Afuresertib, AZD2014, AZD5153, Dactolisib, Erlotinib, Gefitinib, GNE-317, I-BRD9, Ibrutinib, Lapatinib, MK-2206, OF-1, OSI-027, Pictilisib, Uprosertib and Wnt-C59) showed a significant positive correlation with patient risk scores, suggesting that these drugs may not be effective in patients with high-risk ESCA. Conversely, the IC50 values of 3 drugs (Staurosporine, Dasatinib, and Luminespib) exhibited a strong inverse relationship with patient risk scores, suggesting these medications could be beneficial for high-risk ESCA patients. 7