Anoikis classification of lung squamous cell carcinoma reveals correlation with clinical prognosis and immune characteristics

We used publicly available datasets in this study. 504 clinical samples from TCGA, which include 51 normal samples and 502 tumor samples, with one sample having been excluded. 43 LUSC samples were sourced from the GSE50081↗ dataset, 238 samples from GSE30219↗, GSE73403↗ and GSE37745↗. 49976 transcriptome profiles data retrieved from GDC TCGA. The patients diagnosed with LUSC were characterized by various clinical factors such as age, sex, tumor grade, TMN staging, duration of survival, and their status regarding survival. For this analysis, this data is secondary data and does not contain any data which can identify individual. Therefore, ethics statement has been exempted by the Ethical Review Committee of Tianjin Medical University Cancer Institute and Hospital.

We employed the BioManager and limma packages to analyze the differentially expressed genes (DEGs) in normal and tumor samples in anoikis–LUSC TCGA data. The screening criteria of DEGs is |log2(fold change)| >1 and adjusted FDR < 0.05. We selected the top 50 up-regulated and 50 down-regulated genes for our study. The heatmap and volcano plot was performed with ‘limma’ and ‘pheatmap’ packages of R software. We used the survival package, limma package, survminer package, reshape2 package, and RColorBrewer package in R. Adjusted p < 0.05 was considered statistically significant. Anoikis–TCGA and GEO clinical data were obtained from uniCOX analysis expressed by forest, and the data expression through prognostic networks.

A total of 717 genes associated with anoikis were identified and categorized through the official GeneCards website (https://genecards.org/↗) and the Harmonizome platform (https://maayanlab.cloud/Harmonizome/↗). The genomic data for TCGA LUSC were sourced from the official UCSC Xena website (https://xena.ucsc.edu/↗). The investigation analyzed copy number variations (CNV), genetic mapping, and the somatic mutation rates for the identified anoikis-related factors. The mutation frequency derived from CNV was computed utilizing a Perl script algorithm, and the data was visually represented using the ‘RCircos’ tool.

We used the software package ‘Consensus Cluster Plus’ R to perform a consistency analysis and clustered and scored the data according to the expression profiles of anoikis-related genes. We used the k-means algorithm to classify TCGA–LUSC and GSE50081↗ data into several clusters. According to the consensus matrix, consistent cluster score, and cumulative distribution function (CDF) curve, we determined the maximum subtype number k (k = 2) and thoroughly assessed the optimal number of clusters.

We conducted the analysis of genomic variations to inquire the functional relationships between different groups. The data of ‘c2.cp.kegg.symbols.gmt↗ ’ and ‘c5.go.symbol s.gmt↗ ’ from MSigDB could be used to functional enrichment analysis. Genomic pathways could be identified by using the R ‘GSVA’ package. Adjusted p < 0.05 was considered as statistically significant. We used R ‘GSEABase’ and ‘GSVA’ packages to achieve the immune assessment. We used the R ‘limma’ package to analyze immune-related genes between clusters and used the R ‘ggpubr’ and ‘heatmap’ packages to draw the boxplot and heat map analysis results. We used R ‘clusterProfiler’ and ‘enrichplot’ packages for GO/KEGG enrichment analysis and results visualization respectively.

We established a model to assess the association between the LUSC prognosis and the anoikis. Firstly, we applied univariate analysis to recognize survival-related genes in differentially expressed genes (DEGs) based on the cluster. Secondly, an anoikis-associated prognostic risk model was been developed by using Lasso–Cox analysis. Next, we divided all the LUSC samples into risk train set (n = 269) and risk test set (n = 269) randomly, and calculated the anoikis risk score by the model equation, which according to the three essential gene expression profiles (FADD, BAG4, SNAI1). Model Formula: Risk score =FADD* 0.156653387181982-BAG4* 0.194369656312506 + SNAI1* 0.174563712695362. We used Kaplan–Meier survival analysis to evaluate the effect of risk score in clinical prognosis between high-risk and low-risk groups, which were defined according to the standard classification rules. The accuracy of the 1-, 3-, and 5-year predictions for the train, test, and all groups was verified by receiver-operator characteristic (ROC) curves. Principal component analysis (PCA) was performed based on the expression of the anoikis-related characteristic genes. The distribution of the groups was explored using t-distributed stochastic neighbor (t-SNE) analysis using the ‘Rtsne’ ‘R package’. Additionally, the dataset GSE30219↗, GSE73403↗ and GSE37745↗ serves to validate the anoikis model.

We processed the data by constructing independent prognostic analyses using the R ‘survival’ and ‘survminer’ packages. Adjusted p < 0.05 was considered as statistically significant. Forest plots plotted the results of the independent prognostic analyses. The entire risk expression file and typing results were expressed through the heatmap package by drawing the risk heatmap. We constructed a nomogram based on the multiple clinical characteristics and risk scores of LUSC patients, which can show the corresponding 1-, 3- and 5-year survival probabilities based on the different risk scores. We used the ‘rms’ R package to make this nomogram, and used the ‘ggpubr’ R package to draw the Box plot. Observing the relationship between sample typing and risk scores through the Sankey diagram to express the accuracy of patient survival in 1-, 3- and 5-year periods by drawing calibration curves. We used the ‘ggDCA’ R package to construct decision curves in 1-, 3- and 5-year periods.

We used the CIBERSORT algorithm to assess the characteristics of immune cell infiltration within the tumor microenvironment (TME) (https://cibersort.stanford.edu/about.php↗), which categorizes different immune cell types. These findings from the CIBERSORT analysis formed the basis for the subsequent evaluations. We examined the immune cell profiles within the anoikis risk classification. Utilizing the ‘limma’ R package, we assessed the heatmap illustrating the relationships between risk scores and immune checkpoint genes. The data regarding immune cell infiltration and risk scores were processed using the R package ‘limma’. The disparities between high-risk and low-risk groups were illustrated through bar charts created with the ‘ggpubr’ package, correlation visualizations utilizing the ‘ggpubr’ library for ordered gene expression data.

Moreover, we assessed the variations in TME scores and the expression of immune-related genes across the two risk categories. Subsequently, we estimated the tumor-infiltrating immune cell (TIC) abundance profile and immune-related biological functions using a single-sample (ssGSEA) algorithm. Adjusted p < 0.05 was considered as statistically significant.

Thee violin diagrams were used to visually demonstrate the differences between stromalScore, ImmuneScore, and ESTIMATEScore risk groups. We used expression data and sensitivity data of targeted drugs from Genetics of Drug Sensitivity in Cancer (GDSC)(https://www.cancerrxgene.org/↗). To obtain score files of the sensitivity of the train and test groups to the different therapeutic drugs, we firstly analyzed the differential expression data of the anoikis gene in TCGA and GEO for drug sensitivity by using R software including ‘limma’, ‘oncoPredict’, and ‘parallel’ packages. And then, the obtained drug susceptibility score files and sample risk files were processed using the ‘ggplot’ and ‘ggpur’ software packages, and the drug susceptibility results for the high-risk and low-risk groups were plotted draw box plots. Multifactorial results will be available, including FADD, SNAI1, and BAG4.

Independent prognostic genes were identified through the analysis of univariate and multivariate Cox regression, and forest plot were used to display the variate (p value, hazard ratios (HR), and 95% confidence intervals (CI). The Kaplan–Meier survival analysis can enable to show prognostic difference between high and low levels of independent prognostic gene expression. p < 0.05 was considered statistically significant.

The LUSC cell lines HCC95/NCI H520 and 293 T were purchased from the American Type Culture Collection (ATCC). The lentiviral particles were generated using 293 T cell lines. The sequence of sgRNA that targets SNAI1 is as follows: SNAI1-sg2: ctctgaggccaaggatctcc; The alternative sgRNA designed for SNAI1 is noted as SNAI1-sg6: gaagtagaggagaaggacga. Following the knockout of SNAI1, we tested the genome indel, mRNA relative expression, proliferation, and apoptosis in cells.

For this study, we use p values and p < 0.05 was considered to be statistically significant (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001). We used a univariate Cox proportional risk regression model to calculate risk ratios for univariate analyses and a multivariate Cox regression model to determine independent prognostic factors. Significance of differences in risk scores at different stages assessed by Wilcox test. All analyses were conducted with R version 4.3.0 in the present study. Related R packages including ‘survival’, ‘BioManager’, ‘limma’, ‘heatmap’, ‘survminer’, ‘reshape2’, ‘RColorBrewer’, ‘RCIRCOS’, ‘consensus cluster plus’, ‘Rtsne’, ‘umap’, ‘ggplot’, ‘ggpubr’, ‘pheatmap’, ‘GSEABase’, ‘GSVA’, ‘GSEBase’, ‘clusterProfiler’, ‘enrichplot’, ‘caret’, ‘glmnet’, ‘timeROC’, ‘ggplot2’, ‘ggalluvial’, ‘dplyr’, ‘regplot’, ‘rms’, ‘ggDCA’, ‘corrplot’, ‘vioplot’, ‘tidyverse’ and ‘ggExtra’ were downloaded from the Lanzhou University Open-Source Society.

The investigation into the significance of genes associated with anoikis in the context of LUSC was conducted. The study provided insights into genes’ differential expression and subsequent transcriptional activity. A comprehensive analysis revealed 407 genes associated with anoikis, identified from the differences between the normal and diseased cohorts in the TCGA dataset. then we selected 197 differentially expressed genes (DEGs) and used the thermal map shows that 50 ARGs were up-regulated and 50 ARGs were down-regulated genes between anoikis -normal and anoikis-tumor subtypes based on the screening criteria of |logFC(fold change)|>1, FDR < 0.05 (Figure 1A). The volcanic map shows the distribution of these ARGs (Figure 1B; Supplementary Table 1↗). A total of 19 out of 197 differentially expressed genes (DEGs) were identified as related to overall survival (OS) in the univariate Cox regression analysis, as illustrated in the forest plot (Figure 1C; Supplementary Table 2↗). The study of anoikis-related genes by network graph displays the association of risk factors based on p values (Figure 1D). Regarding copy number variations, a higher amplification frequency was observed in 10 genes than in deletions. Specifically, CLDN1, BAG4, FADD, MUC1, PDK4, SERPINA1, CHEK2, HGF, GDF2, and SNAI1 showed extensive CNV amplification. in contrast, several genes like LATS1、ERBB4、SPINK1、TJP3、CD151、and THBS1 were found to display deletions in their CNV (Figure 1E). The graphical representation of changes in CNV across 19 anoikis-related genes is illustrated on the chromosomes (Figure 1F).

From the TCGA–LUSC and GEO-GSE50081↗ datasets, a total of 19 factors related to anoikis were identified. A consensus clustering algorithm was employed to categorize patients with LUSC based on the expression data of the 19 identified genes. To explore the expression patterns of genes associated with anoikis in LUSC, a total of 544 samples derived from 501 cases in TCGA–LUSC along with 43 samples from GEO-GSE50081↗ were analyzed using the consensus clusterplus approach, focusing on the expression data of 19 anoikis-related genes. According to our findings, k = 2 was the best grouping variable for dividing the dataset into Cluster A and Cluster B (Figure 2A). Cluster analysis revealed two distinct groups of anoikis-associated genes, designated as Group A with 220 samples and Group B consisting of 324 samples. Dimensionality reduction techniques such as PCA, tSNE, and UAMP analysis indicated a pronounced clustering of samples from types A and B (Figure 2B–D). According to the survival analysis, the survival results of the two clusters were significantly different (p = 0.005) (Figure 2E).

Various assessments indicated that these four genes exhibited elevated expression levels in cluster B, as supported by the heatmap analysis (Figure 3A,B). Additionally, two distinct subtypes exhibited a notable accumulation of immune cells (Figure 3C). A GSVA analysis was performed to compare the biological activities between clusters A and B (Figure 3D). In cluster B, significant enrichment was observed for pathways such as KEGG_PROGESTERONE_MEDIATED_OOCYTE_MATURATION. KEGG_BASAL_TRANS­CRIPTION_FACTORS, KEGG_CELL_CYCLE, KEGG_HOMO­LOGOUS_RECOMBINATION. KEGG_DNA_REPLICATION. and KEGG_MISMATCH_REPAIR, while the other pathways showed strong enrichment in cluster A (Figure 3D). In addition, an analysis of the pathways involved in clusters A and B was conducted using Gene Set Enrichment Analysis (GSEA), revealing that the genes significantly expressed in cluster A were predominantly associated with KEGG_COMPLEMENT_AND_COAGULATION_CASCADES, KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION. KEGG_HEMATOPOIETIC_CELL_LINEAGE, and KEGG_SYS­TEMIC_LUPUS_ERYTHEMATOSUS (Figure 3E). The analysis indicated that KEGG_METABOLISM_OF_XENOBIOTICS_BY_CYTOCHROME_P450 was notably prevalent in Cluster B (Figure 3F).

The categorization of anoikis holds significant relevance for the clinical outcomes associated with LUSC. A prognostic model was developed as an initial step to gain deeper insights into the predictive traits of anoikis. A total of 717 genes associated with anoikis were discovered from the clusters, and through univariate Cox analysis, 19 genes linked to prognosis were identified. Secondly, Lasso–Cox regression analysis was conducted on the 19 prognostic genes identified (Figure 4A,B), and three significant genes were discovered (FADD, SNAI1, and BAG4 (Supplementary Table 3↗). The risk score prognostic model were constructed according to the three anoikis-related genes. Risk score=FADD* 0.156653387181982-BAG4* 0.194369656312506 + SNAI1* 0.174563712695362. Then, The LUSC samples were classified into high-risk and low-risk categories according to their risk scores, revealing that the prognosis for those in the high anoikis risk category was less favorable than that of the low anoikis risk category (Figure 4C, p< 0.001). This result was confirmed by Kaplan–Meier curves of the train and test groups (Figure 4D–E, p = 0.004, p = 0.015). We also confirmed the high accuracy of the anoikis risk score in predicting prognosis after 1, 3, and 5 years. The area under the curve (AUC) for the total sample was 0.576 over 3 and 5 years. 0.642 and 0.639 for 1, 3, and 5 years (Figure 4F). The AUC in the training cohort was determined to be 0.532, 0.678, and 0.655 for 1, 3, and 5 years (Figure 4G). The AUC values recorded in the evaluation set were 0.617, 0.606, and 0.629 for 1, 3, and 5 years (Figure 4H). Additionally, we used the external cohort (the dataset GSE30219↗, GSE73403↗, and GSE37745↗) to validate the anoikis model further. The results also showed the prognostic ability of the anoikis model and were statistically significant (Supplementary Figure 1↗). Considering the age in the analysis, Gender, T N M classification, the risk score was associated with OS (HR = 1.984, 95% confidence interval 1.356–2.90, p < 0.001) (Figure 4I). The risk heatmap showed that BAG4 was a low-risk gene, and FADD and SNAI1 were identified as genes associated with high risk (Figure 4J). The analysis of differences reveals distinct characteristics between the patients in clusters A and B, highlighting the relationship between patient categories and their outcomes (Figure 4K). Ultimately, the findings highlight the promising ability of the anoikis risk score as a prognostic indicator. The alluvial diagram revealed that the patients in the low-risk category and most of those who survived predominantly belong to subtype B. In contrast, those exhibiting elevated risk scores are classified under subtype A, consistent with earlier findings (Figure 4L).

We integrated the risk scores and various clinical and pathological variables to enhance the prediction accuracy for the LUSC prognosis. Utilizing the developed nomogram, the overall scores derived from the clinical data can yield a survival rate percentage for patients at 1, 3, and 5 years (Figure 5A). The survival rates for patients at the 1, 3, and 5-year marks were determined to be 0.919, 0.783, and 0.704. The calibration graph indicates that the predicted outcomes for the 1-year, 3-year, and 5-year overall survival rates based on the nomogram closely aligned with the ideal scenario for the entire cohort (Figure 5B). The survival curves indicate that individuals classified in the high-risk category faced a heightened risk compared to their low-risk counterparts (Figure 5C). The DCA analysis revealed that the developed nomogram outperformed other clinicopathological variables in forecasting survival outcomes for patients with LUSC (Figure 5D–F).

To investigate the association between the infiltration of innate immune cells (ICs) and anoikis-related genes, we employed the CIBERSORT algorithm to assess the distribution of these immune cells (Figure 6A,B). The analysis presented in the heatmap illustrates how ICs are associated with the anoikis genes (Figure 6C,D). A notable positive association was observed between the risk score and regulatory T cells (Tregs) and M2 macrophages (Figure 6E,F). Analysis of the differential scores within the tumor microenvironment (TME) revealed notable disparities between the groups classified as high-risk and low-risk regarding StromalScore, ImmuneScore, and ESTIMATEScore. Specifically, the TME score was markedly elevated in the high-risk cohort compared to the low-risk cohort (Figure 6G).

The use of pharmacological treatments is vital in managing tumors. Subsequently, an analysis was carried out regarding the sensitivity of drugs and the prognostic biomarkers in patients with LUSC, revealing notable differences in 30 drugs between the high-risk and low-risk categories. The IC50 values for a selection of 15 pharmacological agents (Afatinib, Afuresertib, AZD3759, BI-2536, BMS-345541, Dihydrorotenone, Erlotinib, Gefitinib, Lapatinib, ML323, OSI-027, Pevonedistat, Pyridostatin, SB505124, Uprosertib), were found to be elevated in individuals identified as high-risk, suggesting that the vulnerability to these medications in the high-risk category was less compared to that observed in the low-risk category (Figure 7A), with the IC50 values for the 15 studied medications (AZD8186, BMS-536924, BMS-754807, Dasatinib, Epirubicin, ERK-2440, Foretinib, GSK269962A, IGF1R-3801, JAK-8517, JQ1, Ribociclib, Trametinib, WIKI4, WZ4003) were markedly reduced in individuals classified as high-risk, indicating that these medications have a strong efficacy in individuals categorized as high risk (Figure 7B).

Subsequently, univariate and multivariate Cox regression analysis was conducted on the expressions of the three prognostic genes and clinical factors, focusing on p values, risk coefficients (HR), and confidence intervals (CI). The analysis indicated that BAG4, SNAI1, and the T category of TNM classification may serve as independent prognostic indicators for LUSC (Figure 8A,B). The study of the Kaplan–Meier curves indicated that individuals with high levels of SNAI1 expression experienced a poorer prognosis than those with lower levels, with the findings reaching statistical significance (Figure 8C,D).

The SNAI1 gene was disrupted in the HCC95 and NCI H520 cell line using CRISPR technology with sgRNA to assess the prognostic value of SNAI1 in lung squamous cell carcinoma (LUSC). We tested the genome indel efficiency by sanger sequencing (Figure 9A), and tested the SNAI1 mRNA relative expression by qPCR (Figure 9B). The proliferation was inhibited prominently and apoptosis was significantly increased in cells of knocking out the SNAI1 (Figure 9C,D). In conclusion, these findings indicate that the prognostic gene SNAI1 is crucial for the proliferation of LUSC.