Anoikis-Related Genes Signature Contributes to Predicting Prognosis and Response to Immunotherapy in Lung Squamous Cell Carcinoma

Introduction

According to the 2025 CA Cancer Statistics, over 680 000 lung cancer survivors live in the United States, yet the overall 5-year survival rate remains a mere 27%, with squamous cell carcinoma constituting a major histological subtype [1]. Lung squamous cell carcinoma and lung adenocarcinoma are often grouped as non-small cell lung cancer, but they harbor distinct molecular profiles [2]. LUSC is more aggressive than LUAD, with larger tumors, higher clinical stage, and a higher male-to-female ratio [3]. Lung squamous cell carcinoma has a worse prognosis than adenocarcinoma and lacks effective targeted therapies [4]. Despite the clinical promise of immune-checkpoint inhibitors and targeted therapies, most NSCLC patients fail to achieve durable responses, highlighting the critical role of the tumor micro-environment (TME) and its spatial heterogeneity in shaping therapeutic outcomes. Immune-excluded niches, characterized by fibrosis, driven by cancer-associated fibroblasts, hypoxia, and macrophage-mediated extracellular matrix remodeling, foster immune exclusion and resistance to immunotherapy. This process links epithelial-to-mesenchymal transition (EMT) and anoikis resistance to an immunologically cold tumor micro-environment, presenting a major barrier to treatment [5].

Anoikis is a form of programmed cell death that occurs when anchorage-dependent cells detach from the extracellular matrix or fail to form proper cell-matrix adhesions. This loss of attachment triggers apoptotic pathways, primarily through integrin signaling disruptions. Anoikis plays a crucial role in tissue homeostasis by eliminating detached or misplaced cells. Resistance to anoikis is a key mechanism in cancer metastasis, allowing tumor cells to survive and spread to distant sites. However, cancer cells frequently acquire anoikis resistance, enabling anchorage-independent growth, survival in circulation, and metastatic colonization [6]. The persistence of recurrence, even after complete surgical resection, underscores the role of intrinsic mechanisms like EMT, stemness, and anoikis resistance in determining long-term outcomes [7,8]. EMT is a biological process in which epithelial cells, which are typically characterized by tight cell–cell junctions, apical-basal polarity, and limited migratory ability, undergo a series of molecular changes that allow them to acquire mesenchymal traits. These features include increased motility, loss of cell polarity, and enhanced migratory and invasive capabilities. EMT plays a key role in cancer metastasis. Tumor cells undergoing EMT gain the ability to detach from the primary tumor and to invade surrounding tissues, and facilitates their spread to distant organs. Cancer stem cells (CSCs), characterized by self-renewal capacity and resistance to chemotherapy, play a pivotal role in cancer relapse. Higher frequencies of CSCs in locally advanced NSCLC are associated with an increased risk of relapse [9]. ECM remodeling is a hallmark of tumor progression, enhancing cancer cell migration and anoikis resistance. Matrix metalloproteinases (MMPs), such as MMP-2 and MMP-9, are critical in this process and facilitate metastatic spread. Immune-checkpoint inhibitors (ICIs), including monoclonal antibodies targeting PD-1, PD-L1, and CTLA-4, have revolutionized the treatment of advanced LUSC. However, immune-related adverse events (IRAEs), a unique class of toxicities associated with ICIs, commonly affect the skin, gastrointestinal system, and endocrine organs. The frequency and severity of IRAEs are heightened when ICIs are combined with chemotherapy [10]. Despite the success of ICIs, resistance to these agents remains a significant challenge in the treatment of lung cancer, particularly in LUSC, with many patients experiencing disease progression due to acquired resistance. Emerging bispecific antibodies and immune-cell engagers show promise in overcoming ICI resistance in LUSC [11,12]. The present study explored the complex interplay between the TME and immune evasion in LUSC, emphasizing the need for novel therapeutic strategies that target these mechanisms to improve patient outcomes. It is crucial to understand mechanisms such as stemness, epithelial-to-mesenchymal transition, and anoikis resistance that influence therapeutic responses in LUSC.

Material and Methods

DATA COLLECTION:

This study utilized data from The Cancer Genome Atlas (TCGA) Lung Squamous Cell Carcinoma dataset. Data were downloaded from UCSC Xena [13]. The dataset includes gene expression data (RNA-seq), clinical data (survival data, tumor staging, clinical characteristics), and molecular features of LUSC patients. Gene expression data was pre-processed and normalized to ensure uniformity and comparability. Clinical data of LUSC patients, including age, sex, smoking history, tumor staging, and treatment regimens, were integrated with gene expression data for survival analysis and immune-cell infiltration analysis. Additionally, 2 LUSC-related datasets from the Gene Expression Omnibus (GEO) database – GSE33479 and GSE12472 – were used [14]. All GEO datasets were log2 transformed for standardization to ensure compatibility with TCGA data. EMT, stemness, and anoikis-related genes were obtained from GeneCards. Patient data from TCGA and GEO are publicly available, de-identified data, so no IRB approval or informed consent was required.

DIFFERENTLY EXPRESSED GENES:

Differentially-expressed genes in LUSC were identified using the limma R package. A Venn diagram was used to obtain the common genes in the 3 datasets. The expression of these genes was compared between the normal and tumor samples. A PCA graph was constructed using the Rtsne package. A heatmap was plotted using the pheatmap package. GO, KEGG, GSEA enrichment analysis, and protein interaction network analyses were also performed.

SURVIVAL ANALYSIS:

The association between differentially-expressed gene levels and overall survival in LUSC patients was evaluated using Cox proportional hazards regression and Kaplan-Meier survival curves. Gene expression was categorized into high- and low-expression groups based on the median expression level. The log-rank test was applied to compare the survival curves between high- and low-expression groups. The hazard ratio (HR) for each gene was calculated using Cox regression, along with its 95% confidence interval (CI), to assess the prognostic role of these genes in LUSC survival.

NOMOGRAM CONSTRUCTION:

A prognostic nomogram was constructed using clinical variables (age, smoking history, tumor stage) and differentially-expressed gene data. The nomogram provides a visual tool to predict the 3-year and 5-year survival probabilities of LUSC patients based on their clinical and molecular profiles. Variables were selected based on Cox regression analysis to determine the key prognostic factors in LUSC. These factors were then incorporated into the nomogram model. Each variable’s regression coefficient was used to calculate the total score for each patient, which was converted to the corresponding survival probabilities for 3 and 5 years. Calibration curves for 3-year and 5-year survival were constructed to evaluate the predictive accuracy of the nomogram. The actual survival probabilities were compared with the predicted probabilities. ROC curves were generated to evaluate the performance of the nomogram at predicting survival for 1, 3, and 5 years. The area under the curve (AUC) was calculated to quantify the model’s discrimination ability. Decision curve analysis (DCA) was performed to assess the net benefit of the nomogram at different risk thresholds.

IMMUNE INFILTRATION ANALYSIS:

CIBERSORT is a computational method based on gene expression data to estimate the relative proportions of immune-cell subtypes in tumor tissue. It was applied to TCGA-LUSC dataset to perform immune-cell subtype analysis. The LM22 immune-cell signature matrix was used to estimate the proportions of different immune-cell types (B cells, T cells, macrophages, dendritic cells) in each sample. IOBR packages was utilized [15]. The input data were normalized RNA-seq data. CIBERSORT analysis was conducted to estimate the fraction of different immune-cell types in each tumor sample. In addition to CIBERSORT, quanTiseq and SVR (Support Vector Regression) methods were used to estimate immune-cell infiltration in the tumor samples, providing a quantitative description of immune-cell interactions.

STATISTICAL ANALYSIS:

All statistical analyses were performed using R software (version 4.5.1) and Python (version 3.13.7). Survival analysis was conducted using the survival package, nomogram construction was performed using the rms package, immune-cell infiltration analysis was done using CIBERSORT, and ROC curve analysis was carried out using the pROC package. The log-rank test was used to compare survival curves. For immune-cell infiltration analysis, group differences were evaluated using the Wilcoxon rank-sum test or the Mann-Whitney U test. All statistical tests were 2-sided, and P values <0.05 were considered statistically significant.

Results

Table 1 summarizes the demographic and clinical characteristics of the study population. A total of 321 participants were included across the 3 datasets: 178 from TCGA-LUSC, 108 from GSE33479, and 35 from GSE12472. Most participants were male. The mean ages were similar across datasets, ranging from 62 to 68 years. In TCGA-LUSC, most participants were former smokers, particularly those who had quit within 15 years. Detailed staging information was available for the TCGA and GSE12472 datasets. In TCGA-LUSC, tumor stages I–III were predominant (97 stage I, 38 stage II, 38 stage III), while metastatic disease was rare (3 stage IV). Most tumors were T2 (116 cases) and N0 (117 cases). In GSE12472, tumors were classified into T1–T3, and metastatic patterns included 15 cases without metastasis, 9 with lymph node involvement, and 12 with distant metastases. Histological grading information was available only for the GSE33479 dataset, which included a wide spectrum of normal, premalignant, and malignant epithelial states: 13 normal tissues, 15 hyperplasia, 15 metaplasia, 13 mild dysplasia, 13 moderate dysplasia, 12 severe dysplasia, 13 carcinoma in situ, and 14 invasive squamous cell carcinoma samples.

Figure 1 illustrates a volcano plot comparing the tumor and normal samples from TCGA-LUSC (upper left). The x-axis represents the log2 fold change (log2FC), which quantifies the difference in gene expression between the 2 groups, while the y-axis shows the negative logarithm of the P value, indicating the statistical significance of the differences. P values <0.05 and log2 fold change greater than 2 indicate strong statistical significance. The volcano plot on the upper right compares the expression profiles of tumors versus normal samples from the GSE33479 dataset. The volcano plot on the bottom left displays the results of a comparison between tumors and normal samples from the GSE12472 dataset. On the bottom right, the Venn diagram illustrates the overlap of significant genes across the 3 comparisons: tumor vs normal from the TCGA-LUSC, GSE33479, and GSE12472 dataset. The diagram reveals that 3 genes were significant in all 3 datasets.

In Figure 2, the boxplot on the top left displays the expression levels of S100A7 in lung squamous cell carcinoma tumor versus normal samples. It demonstrates significantly higher expression levels of S100A7 in tumor samples compared to normal tissues (P<0.05). The whiskers represent the range of expression, with outliers indicating extreme values in both tumor and normal samples. The boxplot on the top right illustrates the expression levels of S100A8 in LUSC tumor versus normal samples. The expression was higher in the tumor samples, with a significant difference between the 2 groups (P<0.05). The boxplot on the bottom left shows the expression levels of SPP1 in LUSC tumor versus normal samples. SPP1 is notably upregulated in the tumor samples compared to normal tissues, with a significant difference (P<0.05). The PCA plot provides an overview of the overall gene expression variance in LUSC samples, suggesting a significant difference in the gene expression profiles between these 2 groups.

In Figure 3, the protein–protein interaction network on the top left shows the interactions between the top differentially-expressed genes identified in the datasets. The nodes represent genes, while the edges indicate interactions between them. Colored edges denote different types of interactions such as functional associations, physical interactions, and the size of the nodes reflects the degree of connectivity for each gene. The Gene Ontology dotplot on the top right displays the top 10 enriched GO terms for biological processes, showing the association of DEGs with specific biological functions; regulation of phosphorylation, epidermis development, and cellular response to growth factor stimulus are among the most significant, suggesting that processes related to cellular signaling, proliferation, and development are strongly associated with the differentially-expressed genes. The size of each dot reflects the number of genes associated with the term, and the color intensity indicates the significance. The heatmap on the bottom left illustrates the expression patterns of the top DEGs across the tumor and normal samples. The KEGG barplot on the bottom right displays the top 10 enriched pathways identified by the KEGG database. Pathways such as transcriptional misregulation in cancer, complement and coagulation cascades, and ECM-receptor interaction are among the top enriched pathways.

In Figure 4, the combined GSEA (gene set enrichment analysis) plot shows the running enrichment scores for the top 20 gene sets across the ranked gene list. Each colored line represents a different gene set, with the gene set names and their respective normalized enrichment scores (NES), P values, and false discovery rates (FDR). The KRAS.LUNG_UP.V1_DN gene set has an NES of 1.694, with a P value of 2.46e-03 and an FDR of 2.9e-02. ACT_PATHWAY_UP.V1_DN has an NES of 1.458, a P value of 1.27e-03, and an FDR of 3.8e-02, which enriches the analysis in pathways related to actin dynamics, a critical process for cell motility and cancer metastasis.

In Figure 5, the Kaplan-Meier survival curve on the top shows the overall survival of patients with lung squamous cell carcinoma (LUSC) stratified by S100A7 expression. The hazard ratio (HR) between the high and low groups was 0.96 (95% CI: 0.74–1.23), and the log-rank P value was 7.729e-01. The middle Kaplan-Meier survival curve shows the overall survival of LUSC patients stratified by S100A8 expression. The HR between the 2 groups was 0.99 (95% CI: 0.77–1.27), with a log-rank P value of 9.266e-01. The bottom Kaplan-Meier survival curve presents the overall survival of LUSC patients based on SPP1 expression. The HR was 1.05 (95% CI: 0.82–1.35), and the log-rank P value was 6.881e-01.

In Figure 6, the prognostic nomogram predicts the overall survival (OS) probabilities at multiple time points for patients with lung squamous cell carcinoma. The nomogram integrates multiple clinical and molecular variables, including S100A7, S100A8, SPP1, age, smoking status, and tumor characteristics (T, N, M stage). S100A7, S100A8, and SPP1 expression levels are shown with their corresponding point assignments. Age is shown with a scale ranging from 35 to 85 years. The points in the middle yield a total score for each variable.

In Figure 7, the calibration plot on the top left compared the nomogram predicted probability of 3-year overall survival with the actual 3-year OS observed in the dataset. The calibration plot on the top right illustrates the 5-year OS predictions. The DCA assesses the net benefit of using the nomogram predictions at different risk thresholds for overall survival (OS). The ROC curve evaluates the discriminative ability of the nomogram for predicting OS at 1 year, 3 years, and 5 years. The x-axis represents the false-positive rate, and the y-axis shows the true-positive rate. The AUC values for each time point indicate the performance of the model: AUC at 1 year=0.650, AUC at 3 years=0.632, and AUC at 5 years=0.657. These results indicate that the nomogram has moderate predictive ability, with slight improvement in performance over time.

In Figure 8, the bar plot on the top represents the immune-cell fractions for each sample, as predicted by CIBERSORT. The x-axis shows the proportion of each immune-cell type in each sample. Different colors represent various immune-cell types, such as B cells memory, macrophages M1, neutrophils, T cells CD4 memory activated, and T cells CD8. The color intensity for each bar corresponds to the fraction of that specific immune-cell type in the sample, providing a visualization of immune-cell infiltration levels across all samples. The violin plot on the bottom displays the immune-cell infiltration levels by group for each immune-cell type. The violin shapes indicate the distribution of infiltration scores for each immune-cell type across the high and low groups. There were significant differences in immune-cell infiltration between the high and low groups for T cells, macrophages, tregs and dendritic cells between the 2 groups.

In Figure 9, the bar plot on the top shows the immune-cell fractions for each sample, as predicted by quanTseq. Each bar represents the proportion of different immune-cell types, with each color corresponding to a specific cell type. The relative proportions of each immune-cell type can be seen across different samples, providing insight into the immune composition of each sample. Some samples show high proportions of macrophages M2 and T cells CD4, while others have high proportions of dendritic cells or T cells CD8. The bar plot on the bottom visualizes the immune-cell fractions in a second set of samples, predicted by the SVR method. The immune-cell types are represented by different colors, and the fractions of B cells, T cells, memory-CD4-T, FCGR3A-Mono, naive-CD4-T, dendritic cells, platelet, and NK cells are shown for each sample.

Discussion

We found the immune landscape of LUSC reveals a predominantly immunosuppressive micro-environment, characterized by a high presence of M2 macrophages and tregs, alongside low infiltration of CD8+ T cells. These findings suggest a tumor environment that favors immune evasion, which is typically associated with poor prognosis. The imbalance of a skew towards M2 macrophages, and the significant infiltration of tregs point to a mechanism of immune suppression that can hinder effective anti-tumor responses. Therapies aimed at modulating macrophage polarization or depleting tregs could be beneficial for restoring immune function and enhancing therapeutic efficacy in lung cancer. Comparative analysis using CIBERSORT, quanTseq, and SVR methods showed consistent immune-cell fraction predictions. The correlation between immune-cell infiltration and clinical outcomes highlights the importance of personalized immunotherapy approaches, where immune profiling could guide treatment decisions. Patients with higher CD8+ T cell infiltration may respond better to immune-checkpoint inhibitors, while those with higher M2 macrophage or treg populations might benefit from therapies targeting these suppressive cells. This underscores the potential of personalized immune-based strategies in improving treatment outcomes for LUSC. Our novel approach of integrating anoikis resistance and immune-cell infiltration in predicting the prognosis and therapeutic response in lung squamous cell carcinoma provides a more comprehensive understanding of the tumor micro-environment and its impact on treatment outcomes.

In LUSC, higher immune scores paradoxically correlate with poorer survival, indicating the complexity of the tumor immune micro-environment and the prevalence of immune-suppressive infiltration. Previous research found that transcriptomic profiling identified IGLC7 as a novel immunomodulator in LUSC, correlating with B and T cell activity as well as immune checkpoints [16]. Anoikis, a form of apoptosis triggered by loss of extracellular matrix attachment, acts as a safeguard against metastasis. The acquisition of anoikis resistance, often linked with EMT, enables cancer cells to survive detachment and thrive in circulation, contributing to metastasis and immune evasion [17]. Bit1 suppresses EMT and metastasis by upregulating E-cadherin through antagonism of the TLE1 corepressor, connecting anoikis resistance and EMT regulation in LUSC [7]. In NSCLC, oncogenic pathways, including p130Cas/Src, 14–3–3ζ–mediated suppression of Bim, and loss of Bit1-dependent E-cadherin expression, collectively promote EMT, stemness, and immune evasion [16]. The establishment of radioresistant LUSC cell lines revealed that enhanced DNA repair through activation of the ATM/CHK2 and DNA-PKcs/Ku70 pathways plays a central role in radioresistance, contributing to increased invasiveness and ECM-receptor interaction [18]. USP13 is a transcription factor that drives squamous cell differentiation. It facilitates lineage plasticity by switching club cells to a squamous phenotype, which is critical for LUSC progression. It also stabilizes the MYC protein, enhancing tumor aggressiveness [19]. The activation of squamous differentiation and PI3K/Akt pathways in moderate to severe dysplasias indicates their role in LUSC progression from pre-invasive lesions to invasive cancer [20]. Tumor budding, marked by isolated tumor cells or small clusters, is associated with poor survival, recurrence, lymph node metastasis, and pleural invasion. The link between tumor budding, EMT, and circulating tumor cells plays crucial role in metastasis [21]. Anoikis resistance, which allows cells to survive in low-attachment environments, is a prominent feature of metastasis in LUSC. Bleb formation is a process that enables survival in detachment. Inhibition of bleb formation sensitizes cells to chemotherapy [22]. ITGBL1, overexpressed in LUSC, promotes anoikis resistance and enhances metastatic potential through the AKT/FBLN2 axis [23]. GDH1, regulated by PLAG1, promotes anoikis resistance and tumor metastasis through the CamKK2-AMPK signaling axis in LKB1-deficient lung cancer [24]. mTORC1 activation by substrate stiffness, via the integrin-GSK3β-FTO axis, indicates the role of m6A methylation in anoikis resistance [25]. GRP94 enhances anoikis resistance and peritoneal metastasis via activation of the YAP/TEAD1 pathway [26]. EGFR activation, induced by ADAMTS1-mediated VCAN cleavage, is a key regulator of anoikis resistance and invasion [27]. In LUSC, S100A7 is overexpressed and has been implicated in the promotion of anoikis resistance, which enhances the metastatic potential of LUSC [28]. S100A8 is particularly involved in promoting stemness by interacting with key inflammatory pathways, including NF-κB and MAPK, which are known to enhance the survival and self-renewal capacity of cancer stem cells [29]. SPP1, a multifunctional extracellular matrix protein, mediates anoikis resistance in LUSC by activating PI3K-Akt-mTOR signaling, making it a potential therapeutic target to overcome metastatic progression and chemoresistance [30]. The crosstalk between S100A7, S100A8, and SPP1 shows a crucial mechanism driving the metastatic potential of lung squamous cell carcinoma; they collaborate to establish a micro-environment that supports tumor progression, immune evasion, and metastasis. The interplay enables LUSC cells to maintain their stem-like properties and ability to resist apoptosis, and exhibit enhanced invasive and metastatic behaviors. Targeting these pathways might improve patient outcomes in LUSC. Future prospective trials are needed to validate the prognostic model and nomogram in real-world settings to further strengthen its clinical relevance.

Conclusions

We found a significant impact of anoikis resistance and immune-cell infiltration in shaping the prognosis and therapeutic response in LUSC. The prognostic model effectively predicts clinical outcomes. The findings emphasize the need for targeted therapeutic strategies that modulate the tumor micro-environment, particularly immune-cell infiltration and anoikis resistance pathways. Targeting these pathways could improve clinical outcomes by overcoming the immune evasion mechanisms and enhancing the efficacy of immuno-therapeutic interventions in LUSC.