Establishment and validation of a novel anoikis-related prognostic signature of lung squamous cell carcinoma
Highlight box
Key findings
• We developed and validated a novel six-gene anoikis-related prognostic signature that effectively stratifies overall survival in lung squamous cell carcinoma (LUSC) patients across multiple cohorts.
• We identified CSNK2A1 as a key driver of LUSC progression and its knockdown significantly inhibits tumor cell proliferation, migration, and invasion in vitro.
What is known and what is new?
• Anoikis resistance is a hallmark of cancer metastasis, enabling tumor cells to survive after detaching from the extracellular matrix. However, its specific prognostic landscape and impact on the immune microenvironment in LUSC remain poorly defined.
• This study establishes the first anoikis-specific prognostic model for LUSC. It moves beyond simple survival prediction by linking anoikis-related genes (ARGs) to specific immune evasion phenotypes and identifying potential therapeutic vulnerabilities for high-risk individuals.
What is the implication, and what should change now?
• The ARG-signature serves as a robust molecular tool to complement traditional tumor, node, metastasis (TNM) staging, offering more precise risk stratification for LUSC. CSNK2A1 emerges as a promising therapeutic target for mitigating LUSC metastasis.
• Clinicians should consider incorporating anoikis-related biomarkers into LUSC prognostic assessments. Future clinical trials are warranted to evaluate the efficacy of combining Wee1 inhibitors with standard therapies for patients identified as high-risk by this signature.
Introduction
Lung squamous cell carcinoma (LUSC), a major histological subtype of non-small cell lung cancer (NSCLC), constitutes approximately 20–30% of all lung malignancies globally (1,2). Despite notable advancements in multimodal therapies, including targeted therapies and immunotherapies (3-7), the 5-year overall survival rate for patients with advanced LUSC regrettably remains below 15% (8). This poor prognosis is primarily attributable to significant tumor heterogeneity and the high propensity for metastasis (9,10). Existing prognostic models, such as tumor, node, metastasis (TNM) staging and histopathological grading (11), are often insufficient in capturing the comprehensive molecular complexity that underlies disease progression. This limitation underscores the critical need for novel biomarkers that accurately reflect fundamental biological processes and facilitate the implementation of personalized therapeutic regimens.
Anoikis, a specific form of programmed cell death (PCD) triggered by cell detachment from the extracellular matrix (ECM), constitutes a critical intrinsic barrier against metastatic dissemination (12). Cancer cells that acquire anoikis resistance can survive under anchorage-independent conditions, which significantly facilitates their subsequent dissemination to distant organs. Emerging evidence highlights the context-dependent and dual role of anoikis-related genes (ARGs) in tumorigenesis (13). While some ARGs function as tumor suppressors by promoting cell death upon ECM detachment, others paradoxically enhance metastatic potential via mechanisms such as epithelial-mesenchymal transition (EMT) and the activation of survival pathways (14-16). Specifically in LUSC, dysregulated anoikis signaling has been implicated in chemoresistance and immune evasion (17-19). However, systematic investigations into its prognostic implications remain notably scarce. Moreover, previous studies have largely focused on single ARGs or limited pathway analyses (14,20), thereby overlooking the potential combinatorial effects of ARG networks that are likely to more accurately reflect the underlying tumor biology.
Recent advances in multi-omics profiling and machine learning algorithms have enabled the construction of robust gene signatures that effectively integrate multidimensional molecular data (21). However, existing prognostic models for LUSC rarely incorporate anoikis-related mechanisms, despite their established and critical role in metastasis—the major cause of mortality in this malignancy. Furthermore, the interplay between anoikis pathways and the tumor immune microenvironment (TME) remains poorly characterized in LUSC. This existing knowledge gap is particularly critical given the growing recognition of immune checkpoint inhibitors (ICIs) as frontline therapeutic agents, considering that anoikis-related pathways are hypothesized to influence both therapeutic response and intrinsic resistance mechanisms.
Consequently, we hypothesize that a systematic analysis of ARG expression patterns could yield novel prognostic biomarkers and potential therapeutic targets for LUSC. Leveraging transcriptomic data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) cohorts, we developed and employed a robust computational framework to identify ARGs with significant prognostic value. The resulting ARG signature was subsequently validated for its capability to effectively stratify LUSC patients into distinct risk groups exhibiting differential overall survival, metastatic potential, and tumor immune infiltration profiles.
This work successfully establishes the first anoikis-related prognostic model specifically designed for LUSC, providing a powerful molecular tool to significantly refine patient risk stratification and identify individuals who may best benefit from anoikis-targeted therapeutic intervention. Furthermore, our findings significantly advance the understanding of anoikis biology in LUSC progression and critically lay the groundwork for the development of precision oncology approaches in this highly recalcitrant malignancy. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2476/rc).
Methods
Data acquisition
Transcriptomic data and corresponding clinical information for LUSC patients were retrieved from TCGA and the GEO (GSE3141 and GSE73403). The TCGA-LUSC cohort comprised 496 primary tumor samples and 49 adjacent normal lung tissue samples, with 489 patients possessing complete clinical and survival information. The RNA sequencing (RNA-seq) raw count data and corresponding clinical data for the TCGA-LUSC cohort were downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/). Subsequently, the raw count data were subjected to a Variance Stabilizing Transformation (VST) using the DESeq2 R package (22). Similarly, the microarray data for the GEO datasets (GSE3141, n=52; GSE73403, n=69) were downloaded from the GEO portal (https://www.ncbi.nlm.nih.gov/geo/). These two datasets, containing 52 and 69 LUSC samples, respectively, were processed for batch effect correction using the IOBR package (23) using the ComBat method and subsequently integrated for analysis. The GSE3141 dataset utilized the GPL570 platform [(HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array], whereas the GSE73403 dataset was generated based on the GPL6480 platform (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F).
ARGs were obtained via systematic search across the GeneCards (https://www.genecards.org) and Harmonizome portals (https://www.immport.org) databases. We retained genes with a relevance score greater than 0.4 from the GeneCards database, preserved all genes from the two ARG sets in the Harmonizome database, and merged the genes obtained from these two databases, ultimately yielding a set of 640 ARGs for subsequent analysis. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.
Differential expression analysis
RNA-seq data from TCGA-LUSC tumor and adjacent normal samples were initially preprocessed by filtering out genes with zero expression across all samples. Differentially expressed genes (DEGs) between tumor and normal tissues were subsequently identified using the limma R package (24), employing a linear model that accounted for potential batch effects. Genes with |log2 fold change| >0.5 and adjusted P value <0.05 (Benjamini-Hochberg correction) were considered significant.
Prognostic gene screening
The intersection set of genes between the identified DEGs and the ARGs in the TCGA-LUSC cohort was subjected to univariate Cox proportional hazards regression analysis using the survival R package (25). Genes meeting both prognostic significance threshold [hazard ratio (HR) ≠1, P<0.05] and consistency in risk direction (HR >1 for poor prognosis genes, HR <1 for protective genes) were retained. Prior to survival modeling, Venn analysis was performed using the VennDiagram R package (26) to visually represent the overlap between the DEGs and ARGs.
Consensus clustering analysis of ARGs
Unsupervised consensus clustering was applied to the TCGA-LUSC cohort, based on the expression of the overlapping DEGs and ARGs, using the ConsensusClusterPlus R package (27) to identify intrinsic subgroups. The resulting clusters were validated using dimensionality reduction techniques, specifically t-distributed Stochastic Neighbor Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP). The Rtsne R package (28) and the umap R package were utilized for these respective analyses. Kaplan-Meier survival curves were generated to assess the differential survival outcomes across the identified subgroups and clusters, utilizing the survival and survminer R packages. Differential gene expression profiles among these subgroups were visualized using the pheatmap R package.
Evaluation of immune infiltration
The single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm was employed to characterize the differences in immune cell infiltration and functional pathways among the identified subgroups. Additionally, the ESTIMATE algorithm (29) in the estimate R package was utilized to calculate the Immune Score, Stromal Score, and ESTIMATE Score. To further explore the association between ARGs and immune components in LUSC, the mRNA expression matrix was normalized, and the CIBERSORT deconvolution tool (30) was used to estimate the relative abundance of 22 tumor-infiltrating immune cells (TIICs). Benjamini-Hochberg (false discovery rate) correction was applied for multiple comparisons in immune cell infiltration to ensure the robustness of our findings.
Risk signature construction and validation
The TCGA-LUSC cohort was randomly partitioned into a training set (70% of samples) and an internal testing set (30% of samples) using the caret R package. The merged GEO datasets (GSE3141/GSE73403) were designated as the independent external validation cohort. The prognostic ARGs identified by univariate Cox analysis were subsequently subjected to Least Absolute Shrinkage and Selection Operator (LASSO) regression (31) (with 10-fold cross-validation and α=1) to select the optimal feature set and prevent model overfitting. A multivariate Cox proportional hazards model was then derived using the selected genes to construct the final risk signature. Patients were stratified into high-risk and low-risk groups based on the median risk score as the cutoff value. The prognostic performance and survival differences between the risk groups were rigorously validated through Kaplan-Meier survival analysis and time-dependent receiver operating characteristic (ROC) curve analysis. In the cohort comprising all patients, we employed 1,000-iteration bootstrapping to validate the obtained C-index, which was 0.628, demonstrating that the model maintains stable predictive accuracy and low bias.
Construction and evaluation of the nomogram containing ARGs signature
In the entire LUSC cohort, univariate and subsequent multivariate Cox proportional hazards regression analyses were performed to identify independent prognostic factors. The variables assessed included age, gender, tumor (T) stage, node(N) stage, metastasis (M) stage, and the risk score derived from the ARG signature. A prognostic nomogram incorporating the independent predictors was established using the rms R package. The discriminative ability (assessed via time-dependent ROC curve analysis) and calibration (assessed via calibration curves) of the nomogram were evaluated in the entire LUSC dataset. Subsequently, decision curve analysis (DCA) was performed to quantify the clinical net benefit and usefulness of the nomogram.
Functional enrichment analysis
Gene Set Variation Analysis (GSVA) was performed using the GSVA R package (32) to investigate pathway enrichment differences based on the identified cluster or risk groups. The Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets were obtained from the Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/index.jsp). To explore the underlying biological processes and pathways associated with the ARG-based risk signature, Gene Ontology (GO) function analysis and KEGG pathway analysis were conducted using the clusterProfiler R package. P<0.05 was used as the threshold to define statistically significant enrichment.
Cell culture and transfection
The human LUSC cell lines, SK-1 and H226, were acquired and maintained in RPMI-1640 and DMEM media, respectively, supplemented with 10% fetal bovine serum (FBS) and 1% antibiotic-antimycotic solution. Cells were cultured at 37 ℃ in a humidified incubator with 5% CO2. The small interfering RNAs (siRNAs) targeting CSNK2A1 (si-CSNK2A1) and a non-targeting negative control (si-NC) were commercially purchased. The sequences for the human CSNK2A1-specific siRNAs were 5'-GGTGAGGATAGCCAAGGTT-3' and 5'-GTTTGGATATGTGGAGTTT-3'. Cancer cells were seeded in 6-well plates and transfected at approximately 30% confluency using Polybrene-mediated transfection in serum-free medium. CSNK2A1 siRNAs or the negative control (NC) siRNA were introduced into the respective groups. Knockdown efficiency was subsequently validated by Western blotting and quantitative polymerase chain reaction (qPCR) analysis.
Western blotting
Whole-cell lysates were prepared using RIPA lysis buffer supplemented with a protease and phosphatase inhibitor cocktail. Protein concentration in the supernatant was quantified using the Pierce BCA Protein Assay Kit following centrifugation to remove cellular debris. Protein samples were denatured by boiling, separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE), and then electrophoretically transferred onto polyvinylidene fluoride (PVDF) membranes. Following blocking with skim milk solution, membranes were incubated overnight at 4 ℃ with specific primary antibodies. The membranes were washed three times with Tris-Buffered Saline with Tween 20 (TBST) and subsequently incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies for 1 h at room temperature. β-actin was employed as the internal loading control. To visualize both the target protein and β-actin, the blots were cropped prior to antibody incubation, and stripping and reprobing was utilized when the molecular weights of the target and control proteins were similar.
Migration and invasion assays
The wound healing assay was performed to evaluate the two-dimensional migratory potential of the cancer cells. Transfected cells were seeded in 6-well plates. Upon reaching approximately 80% confluence, a linear scratch was created in the cell monolayer using a sterile 200 µL pipette tip. The scratch area was photographed at 0 and 12 h, and the relative gap closure area was quantified using ImageJ software.
The Transwell invasion assay was conducted using chambers pre-coated with Matrigel to assess the invasive potential. Transfected cells (4×104) were resuspended in 400 µL of serum-free basal medium and placed in the upper chamber. Medium supplemented with 10% FBS was added to the lower chamber as a chemoattractant. After 12 h, non-invasive cells on the upper surface were meticulously removed. Invasive cells on the lower surface were fixed with 4% paraformaldehyde, stained with crystal violet, and subsequently imaged. The number of invasive cells was quantified from five random fields per chamber using ImageJ software.
Cell Counting Kit-8 (CCK-8) assay
Cancer cells were seeded into a 96-well plate at a density of 2,000 cells per well (in 100 µL of culture medium) and incubated overnight for attachment. A baseline measurement (day 0) was established by adding the CCK-8 solution to a control group, followed by incubation in the dark for 1 h, and measuring the absorbance at 450 nm. For the remaining groups, the cells were transfected with the respective siRNAs. Cell viability was measured at 24, 48, 72, and 96 h post-transfection by adding the CCK-8 solution. Following a 1h incubation in the dark, the absorbance [optical density at 450 nm (OD450)] was recorded using a microplate reader.
Colony formation assay
The clonogenic potential was assessed by seeding H226 and SK-1 cells, transfected with the respective siRNAs, into 6-well plates at a low density of 2,000 cells per well. After 14 days of culture, the resulting colonies were fixed with 4% paraformaldehyde and stained with crystal violet. Colony counts were quantified using ImageJ software to evaluate the long-term proliferative capacity of the cells.
Drug sensitivity analysis
Drug sensitivity prediction was performed using the oncoPredict R package (33). The analysis utilized the log2-transformed TCGA-LUSC RNA-seq data in conjunction with the Genomics of Drug Sensitivity in Cancer 2 (GDSC2) pharmacogenomic reference dataset. The prediction model incorporated Empirical Bayes batch correction, filtered to retain the top 20% most variable genes, and generated half maximal inhibitory concentration (IC50) predictions via ridge regression.
Molecular docking
Molecular docking analysis between the Wee1 inhibitor (PubChem CID: 10384072) and the target CSNK2A1 protein (PDB ID: 6QY7, retrieved from the RCSB PDB database) was performed using AutoDock Vina v1.2.3 (34,35). The three-dimensional protein structure was preprocessed using PyMOL to remove co-crystallized ligands and solvent molecules. Ligand structures were downloaded from PubChem and converted to the PDBQT format using Open Babel, along with the assignment of Gasteiger partial charges. A docking grid box was defined, centered on the known adenosine triphosphate (ATP)-binding site of the CSNK2A1 protein, to guide the prediction of binding poses. Docking parameters were set to an exhaustiveness of 20 and 10 independent runs per ligand. The most favorable binding conformation was determined by selecting the pose with the lowest predicted binding energy.
Statistical analysis
All data processing and analysis in this article were conducted utilizing R software version 4.2.1, or GraphPad Prism 10. Continuous variables are presented as mean ± standard deviation, The Wilcoxon rank sum test was used for comparison between the two groups. If not specified, the correlation coefficient between different molecules was calculated by Spearman correlation analysis, and a P value of less than 0.05 was used as the criterion for significant difference.
Results
Clinicopathological characteristics of TCGA-LUSC cohort
The baseline clinicopathological characteristics of the entire patient cohort are summarized in Table 1, and these characteristics, stratified by overall survival status, are presented in Table 2. Comparative analysis of survivors (n=279) and non-survivors (n=210) revealed significant age-related and therapeutic disparities. Deceased patients were older (68.27±8.26 vs. 66.37±8.64 years; P=0.01) and exhibited markedly poorer primary therapy outcomes, with a 19% reduction in complete remission rates (73.6% vs. 92.6%). While TNM staging (pT: P=0.19; pN: P=0.57; overall stage: P=0.07) lacked statistical significance, non-survivors demonstrated trends toward advanced disease, including higher proportions of T3/T4 tumors (22.9% vs. 15.8%) and stage III/IV diagnoses (23.6% vs. 14.8%). Racial disparities approached significance (P=0.056), with Black/African American individuals overrepresented in the deceased cohort (11.2% vs. 4.7%). Smoking status (95.6% vs. 97.1%; P=0.52) and radiation therapy utilization (11.4% vs. 13.3%; P=0.67) showed no survival association, likely reflecting cohort homogeneity in smoking exposure (96.2%) and limited radiotherapy application (12.4%). Collectively, these findings underscore age, primary therapeutic response, and potential sociodemographic factors as critical prognostic determinants, while simultaneously highlighting the limited discriminative power of conventional staging systems in this LUSC population.
Table 1
| Characteristics | Overall (n=489) |
|---|---|
| OS | |
| Alive | 279 (57.1) |
| Dead | 210 (42.9) |
| Age (years) | 67.19 [8.52] |
| Gender | |
| Female | 127 (26.0) |
| Male | 362 (74.0) |
| pM | |
| M0 | 402 (98.3) |
| M1 | 7 (1.7) |
| pN | |
| N0 | 313 (64.7) |
| N1 | 126 (26.0) |
| N2 | 40 (8.3) |
| N3 | 5 (1.0) |
| pT | |
| T1 | 111 (22.7) |
| T2 | 286 (58.5) |
| T3 | 69 (14.1) |
| T4 | 23 (4.7) |
| Stage | |
| Stage I | 239 (49.3) |
| Stage II | 156 (32.2) |
| Stage III | 83 (17.1) |
| Stage IV | 7 (1.4) |
| Radiation therapy | |
| No | 340 (87.9) |
| Yes | 47 (12.1) |
| Primary therapy outcome | |
| Complete remission/response | 269 (85.9) |
| Partial remission/response | 5 (1.6) |
| Progressive disease | 25 (8.0) |
| Stable disease | 14 (4.5) |
| Race | |
| Asian | 9 (2.3) |
| Black or African American | 29 (7.6) |
| White | 346 (90.1) |
| Smoker | |
| No | 18 (3.8) |
| Yes | 459 (96.2) |
Data are presented as mean [SD] or n (%). M, metastasis; N, node; OS, overall survival; pM, pathologic M; pN, pathologic N; pT, pathologic T; SD, standard deviation; T, tumor.
Table 2
| Characteristics | Alive (n=279) | Dead (n=210) | P |
|---|---|---|---|
| Age (years) | 66.37 [8.64] | 68.27 [8.26] | 0.01 |
| Gender | |||
| Female | 78 (28.0) | 49 (23.3) | 0.29 |
| Male | 201 (72.0) | 161 (76.7) | |
| pM | |||
| M0 | 230 (99.1) | 172 (97.2) | 0.25 |
| M1 | 2 (0.9) | 5 (2.8) | |
| pN | |||
| N0 | 184 (66.7) | 129 (62.0) | 0.57 |
| N1 | 70 (25.4) | 56 (26.9) | |
| N2 | 19 (6.9) | 21 (10.1) | |
| N3 | 3 (1.1) | 2 (1.0) | |
| pT | |||
| T1 | 69 (24.7) | 42 (20.0) | 0.19 |
| T2 | 166 (59.5) | 120 (57.1) | |
| T3 | 34 (12.2) | 35 (16.7) | |
| T4 | 10 (3.6) | 13 (6.2) | |
| Stage | |||
| Stage I | 142 (51.3) | 97 (46.6) | 0.07 |
| Stage II | 94 (33.9) | 62 (29.8) | |
| Stage III | 39 (14.1) | 44 (21.2) | |
| Stage IV | 2 (0.7) | 5 (2.4) | |
| Radiation therapy | |||
| No | 203 (88.6) | 137 (86.7) | 0.67 |
| Yes | 26 (11.4) | 21 (13.3) | |
| Primary therapy outcome | |||
| Complete remission/response | 188 (92.6) | 81 (73.6) | <0.001 |
| Partial remission/response | 1 (0.5) | 4 (3.6) | |
| Progressive disease | 7 (3.4) | 18 (16.4) | |
| Stable disease | 7 (3.4) | 7 (6.4) | |
| Race | |||
| Asian | 5 (2.3) | 4 (2.4) | 0.056 |
| Black or African American | 10 (4.7) | 19 (11.2) | |
| White | 199 (93.0) | 147 (86.5) | |
| Smoker | |||
| No | 12 (4.4) | 6 (2.9) | 0.52 |
| Yes | 258 (95.6) | 201 (97.1) |
Data are presented as mean [SD] or n (%). M, metastasis; N, node; pM, pathologic M; pN, pathologic N; pT, pathologic T; SD, standard deviation; T, tumor.
Identification of differentially expressed overlapping ARGs
To identify the subset of ARGs that were differentially expressed, we initially performed differential expression analysis (DEA) on the TCGA-LUSC cohort. A total of 9,745 genes were identified as differentially expressed (Figure 1A), with the top 100 most significant DEGs presented in Figure 1B. Subsequently, we compiled a list of 509 ARGs from the GeneCards database and screened 307 differentially expressed ARGs by overlapping them with the identified DEGs (Figure 1C). Univariate Cox proportional hazards regression analysis was then applied to this subset, ultimately identifying 37 ARGs that were significantly associated with overall survival (Figure 1D).
Identification of anoikis-related subgroups
To identify intrinsic, anoikis-related molecular subgroups, we performed unsupervised consensus clustering on the TCGA-LUSC cohort, utilizing the expression matrix of the overlapping ARGs. The optimal parameter for clustering was determined to be k=2, which yielded two distinct subgroups: Cluster A and Cluster B (Figure 2A). Dimensionality reduction techniques, specifically t-distributed Stochastic Neighbor Embedding (t-SNE) and UMAP, further confirmed the clear separation between Cluster A and Cluster B (Figure 2B,2C). Furthermore, Kaplan-Meier survival analysis demonstrated that Cluster B exhibited a significantly poorer overall survival probability compared to Cluster A (Figure 2D). The heatmap visualization of the overlapping ARGs expression confirmed clear differential expression patterns between the two subgroups (Figure 2E). Additionally, GSVA revealed that the two clusters were significantly enriched in completely distinct biological pathways (Figure 2F). These findings collectively validate that the TCGA-LUSC cohort can be robustly classified into prognostically distinct subgroups based on ARG expression profiles.
Immune microenvironment analysis
We applied the ssGSEA algorithm to explore the level of immune infiltration between the two subgroups. The results indicated that Cluster A had a distinct pattern of immune infiltration compared to Cluster B. Furthermore, Cluster A had a significantly higher abundance of CD56bright natural killer (NK) cell than Cluster B, while Cluster B had a significant abundance of immune cells-including activated B cell, activated CD4 T cell, activated CD8 T cell, activated dendritic cell, eosinophil, gamma delta T cell, immature B cell, immature dendritic cell, myeloid-derived suppressor cell (MDSC), macrophage, mast cell, monocyte, NK T cell, NK cell, neutrophil, plasmacytoid dendritic cell, regulatory T cell, T follicular helper cell, type 1 T helper cell, type 17 T helper cell and type 2 T helper cell than Cluster A (Figure 3A).
Furthermore, the CIBERSORT deconvolution algorithm was utilized to estimate the relative abundances of 22 TIIC subtypes (Figure 3B). Correlation analysis of TIICs revealed notable relationships, such as a significant negative correlation between activated NK cells and resting NK cells (r=−0.49), and a strong positive correlation between T cells CD4 memory activated and T cells CD8 (r=0.63) (Figure 3C). Comparative analysis between subgroups demonstrated that Cluster A exhibited significantly higher proportions of memory B cells, T cells CD4 memory resting, NK cells resting, monocytes, M2 macrophages, mast cells resting, and neutrophils compared to Cluster B. Conversely, Cluster B showed elevated levels of T cells follicular helper and NK cells activated relative to Cluster A (Figure 3D).
The ESTIMATE algorithm further indicated that Cluster B had significantly higher StromalScore, ImmuneScore, and ESTIMATEScore than Cluster A, reflecting a more active immune microenvironment in Cluster B (Figure 3E). Additionally, the Tumor Inflammation Signature (TIS) analysis revealed that Cluster B also displayed a markedly higher TIS score compared to Cluster A (Figure 3F).
These integrated results unequivocally demonstrate that the two ARG-defined subgroups possess significantly distinct immune infiltration characteristics.
Construction of an anoikis-related prognostic risk model
The TCGA-LUSC cohort was randomly divided into training and validation sets at a ratio of 7:3. LASSO regression was applied to reduce dimensionality in the standardized training dataset (Figure 4A). The cross-validation plot identified two optimal lambda values (Figure 4B): the first dashed line (lambda =0.0418) corresponded to the minimum mean cross-validated error, while the second (lambda =0.1060) represented the largest lambda value within one standard error of the minimum. Based on prior literature, the first lambda value was selected. Post-dimensionality reduction, nine key genes were identified (AKT2, CLDN1, CSNK2A1, ERBB4, CHEK2, CD151, TDGF1, BAG4, SNAI1). Subsequent stepwise Cox regression analysis refined the model to six genes (AKT2, CSNK2A1, CD151, TDGF1, BAG4, SNAI1). The risk score was calculated as follows:
Risk score = (AKT2 expression × 0.1649) + (CSNK2A1 expression × 0.2070) + (CD151 expression × 0.1527) + (TDGF1 expression × 0.1596) + (SNAI1 expression × 0.2565) − (BAG4 expression × 0.1372).
Patients in the training set were stratified into high- and low-risk groups based on the median risk score. Differential expression patterns of model genes were observed between subgroups, with Cluster B predominantly enriched in the high-risk group (Figure 4C). AKT2, CD151, TDGF1, and SNAI1 exhibited significantly higher expression in the high-risk group, whereas BAG4 showed lower expression. No significant difference was detected for CSNK2A1 between subgroups (Figure 4D). Elevated risk scores correlated with poorer survival outcomes (Figure 4E). Significant survival disparities between risk groups were validated across training, internal validation, and external validation cohorts (Figure 4F-4H). Time-dependent ROC curves demonstrated AUC values of 0.580, 0.721, and 0.725 for 1-, 3-, and 5-year survival in the training set; 0.606, 0.615, and 0.642 in the internal validation set; and 0.590, 0.657, and 0.728 in the external validation set (Figure 4I-4K).
Development of a clinical nomogram integrating risk score
Univariate and multivariate Cox regression analyses identified clinical T-stage and the ARG-based risk score as the two independent prognostic factors for LUSC patients (Figure 5A,5B). A nomogram incorporating clinical variables and risk score was constructed to predict 1-, 3-, and 5-year survival probabilities (Figure 5C). Calibration curves indicated strong concordance between predicted and observed outcomes (Figure 5D). Risk stratification based on the nomogram highlighted significant divergence in cumulative hazard (P<0.001), with the high-risk group exhibiting rapid escalation in mortality over time, contrasting with the stable trajectory of the low-risk group (Figure 5E). DCA confirmed superior net clinical benefit of the nomogram over standalone clinical parameters across 1-, 3-, and 5-year intervals (Figure 5F-5H).
CSNK2A1 knockdown suppresses malignant phenotypes in LUSC
CSNK2A1, the only gene within the risk signature that lacked statistical significance in differential expression (Figure 4D) and prognostic value (Figure S1A-S1F) in the training set, was selected for further functional investigation in human LUSC cell lines. As a serine/threonine protein kinase, CSNK2A1 represents a highly ‘druggable’ target compared to transcription factors. Re-analysis of the TCGA-LUSC data revealed significantly elevated CSNK2A1 expression in tumor tissues compared to adjacent normal tissues (Figure S1G). Western blot (Figure 6A-6C) and qPCR analyses (Figure 6D,6E) confirmed that CSNK2A1 expression was robustly knocked down in both the SK-1 and H226 cell lines, with both mRNA and protein levels significantly reduced relative to the negative control group. The original image correspondivng to the Western blot and other replicate blots can be found in Figure S2. Transwell assays demonstrated a significant reduction in the migratory and invasive capabilities of both cell lines following CSNK2A1 knockdown (Figure 6F-6J). This suppression of migration was further confirmed by the wound healing assays (Figure 6K-6N). Additionally, cell viability assessed by CCK-8 assays (Figure 6O,6P) and long-term proliferative capacity measured by the colony formation assay (Figure 6Q,6R) were significantly diminished in both cell lines after CSNK2A1 knockdown.
Drug sensitivity profiling and molecular docking
The OncoPredict analysis identified ten pharmacological compounds with statistically significant differences in predicted drug sensitivity between the high- and low-risk groups (Figure 7A-7J). High-risk patients showed sensitivity to Wee1 inhibitor, doramapimod, AZD6738, ibrutinib, gefitinib, Eg5_9814, BI-2536, GSK591, and vinblastine, whereas low-risk patients were potentially sensitive to doramapimod. Cluster-specific drug sensitivity analysis revealed preferential efficacy of Wee1 inhibitor for Cluster B and other agents for Cluster A (Figure S3). Intersection analysis highlighted Wee1 inhibitor as the sole compound effective across both subgroups. Molecular docking analysis was performed to investigate the physical interaction between Wee1 inhibitor and the CSNK2A1 protein. The results confirmed stable binding of the Wee1 Inhibitor within the CSNK2A1 ATP-binding pocket (Figure 7K,7L).
Consensus clustering and risk model interplay
An integrated survival analysis (Figure S4A) demonstrated a synergistic prognostic effect: patients classified as Cluster B/high-risk exhibited the poorest overall survival, while the Cluster A/low-risk group had the most favorable outcomes, with intermediate survival observed for the mixed subgroups. Consistently, Cluster B was significantly associated with higher calculated risk scores and a higher corresponding mortality incidence (Figure S4B,S4C).
Discussion
LUSC, which accounts for 20–30% of NSCLC cases (36), presents a formidable therapeutic challenge due to its intrinsic genomic heterogeneity, limited actionable mutations, and high propensity for early metastasis (37). Unlike lung adenocarcinoma, LUSC currently lacks broad-spectrum targeted therapeutic options, relying primarily on chemotherapy and immunotherapy (36). The dismal 5-year overall survival rate for advanced disease (<15%) (8), underscores the urgent necessity for robust prognostic biomarkers reflecting underlying molecular drivers.
This study reports the first comprehensive integration of an ARG signature into prognostic modeling for LUSC, thereby addressing a critical gap in molecular tools available for predicting outcomes in this recalcitrant malignancy. Our findings reveal that dysregulated anoikis signaling not only drives metastatic potential but also profoundly remodels the tumor microenvironment (TME), offering dual utility for both risk stratification and therapeutic targeting. The six-gene prognostic signature (AKT2, CSNK2A1, CD151, TDGF1, BAG4, SNAI1) demonstrated robust and reliable predictive capacity across multiple independent cohorts, significantly outperforming conventional TNM staging in survival discrimination. Our comparative analysis demonstrated that the ARG-signature achieves superior or comparable performance to existing benchmarks, such as the immune-related signature by Liu et al. (38), particularly in long-term (3-year) survival prediction. This suggests that anoikis-resistance pathways, represented by hub genes like AKT2 and SNAI1 in our model, may more accurately capture the inherent aggressive potential of LUSC cells. Furthermore, DCA confirmed that the ARG-signature provides non-redundant prognostic information, significantly enhancing the net clinical benefit when integrated with traditional TNM staging. This reinforces the role of anoikis-related biomarkers as essential complements to the current molecular landscape of LUSC.
The functional dominance of AKT2 and SNAI1 in our model underscores the centrality of PI3K/AKT signaling and EMT in anoikis evasion. This is highly consistent with their established roles in chemoresistance and the survival of circulating tumor cells (CTCs) (39-41). BAG4 emerged as a protective factor, potentially through its dual regulation of apoptosis and NF-κB-mediated inflammation. This functional dichotomy may account for its context-dependent prognostic roles observed across various cancer types (42-44). The lack of CSNK2A1 expression differences between risk groups, contrasted with its essential role in maintaining malignant phenotypes in vitro, highlights the limitations of transcript-level analyses in capturing post-translational kinase activity. Our experimental validation of CSNK2A1’s pro-metastatic function aligns with recent findings in pancreatic cancer (45) and in lung adenocarcinoma (46), suggesting conserved oncogenic roles across epithelial malignancies.
The divergent immune landscapes identified between Cluster A and Cluster B provide critical insights into potential mechanisms of therapeutic resistance. Cluster B’s combination of elevated stromal and immune scores yet poorer survival may reflect an “inflamed desert” phenotype, characterized by abundant but non-functional immune infiltrates and the dominance of immunosuppressive myeloid populations (47,48). The marked enrichment of exhausted CD8+ T cells and M2 macrophages in this subgroup parallels the immunosuppressive TME patterns observed in ICI-resistant gastric cancer (49). Notably, the high-risk group’s predicted sensitivity to Wee1 inhibitors suggests synthetic lethality opportunities in genomically unstable, anoikis-resistant tumors. The molecular docking validation of Wee1 inhibitor binding to CSNK2A1 introduces intriguing possibilities for dual kinase targeting, though further mechanistic studies are warranted.
Clinically, the nomogram’s superior predictive performance over standalone staging parameters (AUC 0.728 at 5 years) addresses the urgent need for personalized prognostic tools in LUSC. The integration of risk scores with T-stage demonstrates particular utility in stratifying early-stage patients, where current overtreatment/undertreatment rates remain problematic. Our drug sensitivity profiling provides actionable insights for risk-adapted therapy, with high-risk patients potentially benefiting from Wee1 inhibitor/ICI combinations to simultaneously target anoikis resistance and immune evasion.
Despite the promising therapeutic candidates identified, several limitations regarding our pharmacological analysis must be acknowledged. First, the drug sensitivity predictions were performed using the oncoPredict algorithm, which fundamentally relies on transcriptomic data from 2D-cultured cancer cell lines. These models inherently lack the structural and biological complexity of human tumors, such as three-dimensional cell-cell interactions and the influence of the TME. Specifically, the TME in LUSC is characterized by dense stroma and immune infiltration, which can significantly impede drug delivery and modulate intrinsic drug resistance—factors that cannot be fully captured by current computational inferences. Furthermore, while molecular docking provided a structural basis for the interaction between CSNK2A1 and Wee1 inhibitors, these findings remain predictive. Extensive in vivo studies and prospective clinical trials are essential to validate the safety and efficacy of these compounds in LUSC patients stratified by our ARG-signature.
While our multi-cohort validation strengthens generalizability, the retrospective design and absence of treatment-response annotations limit clinical translation. Prospective clinical cohorts are required to confirm the model’s stability and clinical utility in diverse populations. Although we characterized the immune microenvironment and observed high immune exhaustion in the high-risk group, we lacked access to LUSC cohorts specifically treated with ICIs. Consequently, the direct predictive power of our signature for immunotherapy response remains to be validated in future clinical trials. The concentrated experimental focus on CSNK2A1, though mechanistically informative, does not fully capture the signature’s collective biological impact. Future research should employ spatial transcriptomics to resolve the spatial heterogeneity of anoikis signaling and incorporate patient-derived organoid models for robust therapeutic candidate validation. Prospective validation in ICI-treated cohorts will be crucial given the observed TME-modulating effects.
Conclusions
This study definitively establishes anoikis dysregulation as a cornerstone of LUSC progression, providing a robust molecular framework that bridges metastatic biology with immune microenvironment remodeling. The resulting six-gene prognostic model and its associated therapeutic insights offer tangible and actionable advances toward precision oncology in this lethal malignancy.
Acknowledgments
We would like to thank who have devoted a lot to this study, including statisticians, reviewers and editors.
Footnote
Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2476/rc
Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2476/dss
Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2476/prf
Funding: This work was supported by the funds from
Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2476/coif). The authors have no conflicts of interest to declare.
Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.
Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.
References
- Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
- Barta JA, Powell CA, Wisnivesky JP. Global Epidemiology of Lung Cancer. Ann Glob Health 2019;85:8. [Crossref] [PubMed]
- Felip E, Altorki N, Zhou C, et al. Adjuvant atezolizumab after adjuvant chemotherapy in resected stage IB-IIIA non-small-cell lung cancer (IMpower010): a randomised, multicentre, open-label, phase 3 trial. Lancet 2021;398:1344-57. [Crossref] [PubMed]
- Ferone G, Lee MC, Sage J, et al. Cells of origin of lung cancers: lessons from mouse studies. Genes Dev 2020;34:1017-32. [Crossref] [PubMed]
- Antonia SJ, Villegas A, Daniel D, et al. Durvalumab after Chemoradiotherapy in Stage III Non-Small-Cell Lung Cancer. N Engl J Med 2017;377:1919-29. [Crossref] [PubMed]
- Reck M, Rodríguez-Abreu D, Robinson AG, et al. Pembrolizumab versus Chemotherapy for PD-L1-Positive Non-Small-Cell Lung Cancer. N Engl J Med 2016;375:1823-33. [Crossref] [PubMed]
- Rittmeyer A, Barlesi F, Waterkamp D, et al. Atezolizumab versus docetaxel in patients with previously treated non-small-cell lung cancer (OAK): a phase 3, open-label, multicentre randomised controlled trial. Lancet 2017;389:255-65. [Crossref] [PubMed]
- Brahmer J, Reckamp KL, Baas P, et al. Nivolumab versus Docetaxel in Advanced Squamous-Cell Non-Small-Cell Lung Cancer. N Engl J Med 2015;373:123-35. [Crossref] [PubMed]
- Wang C, Yu Q, Song T, et al. The heterogeneous immune landscape between lung adenocarcinoma and squamous carcinoma revealed by single-cell RNA sequencing. Signal Transduct Target Ther 2022;7:289. [Crossref] [PubMed]
- Zhang L, Zhang Y, Wang C, et al. Integrated single-cell RNA sequencing analysis reveals distinct cellular and transcriptional modules associated with survival in lung cancer. Signal Transduct Target Ther 2022;7:9. [Crossref] [PubMed]
- Hendriks LE, Kerr KM, Menis J, et al. Oncogene-addicted metastatic non-small-cell lung cancer: ESMO Clinical Practice Guideline for diagnosis, treatment and follow-up. Ann Oncol 2023;34:339-57. [Crossref] [PubMed]
- Wang Y, Cheng S, Fleishman JS, et al. Targeting anoikis resistance as a strategy for cancer therapy. Drug Resist Updat 2024;75:101099. [Crossref] [PubMed]
- Dai Y, Zhang X, Ou Y, et al. Anoikis resistance--protagonists of breast cancer cells survive and metastasize after ECM detachment. Cell Commun Signal 2023;21:190. [Crossref] [PubMed]
- Dong B, Gu Y, Sun X, et al. Targeting TUBB3 Suppresses Anoikis Resistance and Bone Metastasis in Prostate Cancer. Adv Healthc Mater 2024;13:e2400673. [Crossref] [PubMed]
- Wang X, Gao L, Li H, et al. Integrative analysis of multi-omics data identified PLG as key gene related to Anoikis resistance and immune phenotypes in hepatocellular carcinoma. J Transl Med 2024;22:1104. [Crossref] [PubMed]
- Sun L, Chen Y, Xia L, et al. TRIM69 suppressed the anoikis resistance and metastasis of gastric cancer through ubiquitin‒proteasome-mediated degradation of PRKCD. Oncogene 2023;42:3619-32. [Crossref] [PubMed]
- Primeaux M, Liu X, Gowrikumar S, et al. Claudin-1 interacts with EPHA2 to promote cancer stemness and chemoresistance in colorectal cancer. Cancer Lett 2023;579:216479. [Crossref] [PubMed]
- Chaojun L, Pengping L, Yanjun L, et al. TJP3 promotes T cell immunity escape and chemoresistance in breast cancer: a comprehensive analysis of anoikis-based prognosis prediction and drug sensitivity stratification. Aging (Albany NY) 2023;15:12890-906. [Crossref] [PubMed]
- Foley JM, Scholten DJ 2nd, Monks NR, et al. Anoikis-resistant subpopulations of human osteosarcoma display significant chemoresistance and are sensitive to targeted epigenetic therapies predicted by expression profiling. J Transl Med 2015;13:110. [Crossref] [PubMed]
- Jin L, Chun J, Pan C, et al. The PLAG1-GDH1 Axis Promotes Anoikis Resistance and Tumor Metastasis through CamKK2-AMPK Signaling in LKB1-Deficient Lung Cancer. Mol Cell 2018;69:87-99.e7. [Crossref] [PubMed]
- Swanson K, Wu E, Zhang A, et al. From patterns to patients: Advances in clinical machine learning for cancer diagnosis, prognosis, and treatment. Cell 2023;186:1772-91. [Crossref] [PubMed]
- Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
- Zeng D, Ye Z, Shen R, et al. IOBR: Multi-Omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and Signatures. Front Immunol 2021;12:687975. [Crossref] [PubMed]
- Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
- Farewell V. Modeling Survival Data: Extending the Cox Model. Biometrics 2001;57:652.
- Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 2011;12:35. [Crossref] [PubMed]
- Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-3. [Crossref] [PubMed]
- van der Maaten L. Accelerating t-SNE using Tree-Based Algorithms. Journal of Machine Learning Research 2014;15:3221-45.
- Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
- Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
- Ranstam J, Cook JA. LASSO regression. Br J Surg 2018;105:1348.
- Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
- Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
- Eberhardt J, Santos-Martins D, Tillack AF, et al. AutoDock Vina 1.2.0: New Docking Methods, Expanded Force Field, and Python Bindings. J Chem Inf Model 2021;61:3891-8. [Crossref] [PubMed]
- Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 2010;31:455-61. [Crossref] [PubMed]
- Goldstraw P, Ball D, Jett JR, et al. Non-small-cell lung cancer. Lancet 2011;378:1727-40. [Crossref] [PubMed]
- Lau SCM, Pan Y, Velcheti V, et al. Squamous cell lung cancer: Current landscape and future therapeutic options. Cancer Cell 2022;40:1279-93. [Crossref] [PubMed]
- Liu Z, Wan Y, Qiu Y, et al. Development and validation of a novel immune-related prognostic model in lung squamous cell carcinoma. Int J Med Sci 2020;17:1393-405. [Crossref] [PubMed]
- Cui Y, Lin J, Zuo J, et al. AKT2-knockdown suppressed viability with enhanced apoptosis, and attenuated chemoresistance to temozolomide of human glioblastoma cells in vitro and in vivo. Onco Targets Ther 2015;8:1681-90. [Crossref] [PubMed]
- Fuertes G, Del Valle-Pérez B, Pastor J, et al. Noncanonical Wnt signaling promotes colon tumor growth, chemoresistance and tumor fibroblast activation. EMBO Rep 2023;24:e54895. [Crossref] [PubMed]
- Zmetakova I, Kalinkova L, Smolkova B, et al. A disintegrin and metalloprotease 23 hypermethylation predicts decreased disease-free survival in low-risk breast cancer patients. Cancer Sci 2019;110:1695-704. [Crossref] [PubMed]
- Jiang L, Chen Y, Min G, et al. Bcl2-associated athanogene 4 promotes the invasion and metastasis of gastric cancer cells by activating the PI3K/AKT/NF-κB/ZEB1 axis. Cancer Lett 2021;520:409-21. [Crossref] [PubMed]
- Caprini E, Cristofoletti C, Arcelli D, et al. Identification of key regions and genes important in the pathogenesis of sezary syndrome by combining genomic and expression microarrays. Cancer Res 2009;69:8438-46. [Crossref] [PubMed]
- Bao F, An S, Yang Y, et al. SODD Promotes Lung Cancer Tumorigenesis by Activating the PDK1/AKT and RAF/MEK/ERK Signaling. Genes (Basel) 2023;14:829. [Crossref] [PubMed]
- Zhang Y, Kong R, Yang W, et al. DUSP2 recruits CSNK2A1 to suppress AKT1-mediated apoptosis resistance under hypoxic microenvironment in pancreatic cancer. Cancer Lett 2023;568:216288. [Crossref] [PubMed]
- Wang H, Lv Q, Xu Y, et al. An integrative pharmacogenomics analysis identifies therapeutic targets in KRAS-mutant lung cancer. EBioMedicine 2019;49:106-17. [Crossref] [PubMed]
- Park S, Ock CY, Kim H, et al. Artificial Intelligence-Powered Spatial Analysis of Tumor-Infiltrating Lymphocytes as Complementary Biomarker for Immune Checkpoint Inhibition in Non-Small-Cell Lung Cancer. J Clin Oncol 2022;40:1916-28. [Crossref] [PubMed]
- Salcher S, Sturm G, Horvath L, et al. High-resolution single-cell atlas reveals diversity and plasticity of tissue-resident neutrophils in non-small cell lung cancer. Cancer Cell 2022;40:1503-1520.e8. [Crossref] [PubMed]
- Tsutsumi C, Ohuchida K, Katayama N, et al. Tumor-infiltrating monocytic myeloid-derived suppressor cells contribute to the development of an immunosuppressive tumor microenvironment in gastric cancer. Gastric Cancer 2024;27:248-62. [Crossref] [PubMed]

