We obtained gene expression profiling datasets for two nasal polyp samples associated with chronic rhinosinusitis with nasal polyposis (CRSwNP), specifically GSE136825 (platform GPL20301) and GSE179265 (platform GPL24676), from the NCBI Gene Expression Omnibus (GEO) database. The former dataset was designated as the training set, consisting of 33 control samples and 42 CRSwNP samples, while the latter was used as the validation set, comprising 7 control samples and 27 CRSwNP samples. The raw data from GEO were normalized using the R package “NormalizeBetweenArray.” A total of 338 differentially expressed genes (DEGs) were identified from both published literature and the GeneCards database.

The R package limma was used to identify differentially expressed genes (DEGs) between normal and CRSwNP samples, with a significance threshold of p < 0.05. DEGs with |log2 fold change| > 1 were selected.

We employed the CIBERSORT algorithm to assess the relative proportions of 22 infiltrating immune cell types based on gene expression, with a significance threshold set at p < 0.05. The resulting data was then utilized for further analysis. Group comparisons of immune cell proportions were conducted utilizing the Wilcoxon test. Data visualization, such as histograms, heat maps, and box plots, was created utilizing the “ggplot2” and “vioplot” R packages. Pearson correlation coefficients for each immune cell were calculated using the “corrplot” R package, and the outcomes were illustrated in an accompanying heat map.

We computed the correlation coefficient between the expression of ARGs and the percentage of infiltrating immune cells, and visualized the results using the R package “ggplot.” A significance level of p < 0.05 was applied to detect any significant correlations.

We performed unsupervised cluster analysis on ARGs using the “ConsensusClusterPlus” R package to identify distinct anoikis patterns in CRSwNP. The optimal number of subtypes (k) was determined based on the tendency, smoothness of the cumulative distribution function (CDF) curve, consensus score, and consensus matrix. Principal component analysis (PCA) was conducted using the “ggplot2” R package.

We acquired the “c5.go.symbols” and “c2.cp.kegg.symbols” datasets from the Gene Set Variation Analysis (GSVA) MSigDB database. Subsequently, we employed the R packages “GSVA” and “limma” to assess the modified pathways and biological functions within different clusters linked to ARGs.

We applied the Least Absolute Shrinkage and Selection Operator (LASSO), Support Vector Machine-Recursive Feature Elimination (SVM-RFE), and Random Forest (RF) algorithms to identify key ARGs among the Differentially Expressed Genes (DEGs) from both normal and CRSwNP samples. These selected genes were then compared using a Venn diagram created with the “VennDiagram” R package. Subsequently, ROC curves were constructed to assess the predictive capacity of these signature genes in the training set, and the Area Under the Curve (AUC) was calculated using the “pROC” and “InpROC” R packages, respectively. The predictive efficacy of these signature genes was further validated in the verification set. Finally, a nomogram based on these signature genes was developed using the “rms” R package.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted on the identified genes using the R package “clusterProfiler” to explore the distinct signaling pathways and potential functions associated with the signature genes. Statistical significance was set at p < 0.05.

In the study, we first calculated the correlation coefficients between the expression levels of ARGs and immune infiltrating cells. Subsequently, Spearman’s rank correlation analysis was employed to investigate the association between immune infiltrating cells and the specific genes. Ultimately, Lollipop plots were generated using the R package “ggplot.”

Regulatory networks for miRNA diagnostic biomarkers and transcription factor (TF) diagnostic biomarkers were built utilizing characteristic genes with NetworkAnalyst (http://www.networkanalyst.ca). Furthermore, a file containing interactions between drugs and genes was acquired from the drug-gene interaction database (DGIdb) (https://dgidb.genome.wustl.edu/) and imported into Cytoscape software for visualization.

In order to validate the expression of candidate genes in nasal polyps and inferior turbinate tissues (control groups), total RNA was extracted from the nasal polyps and control groups (three samples for each group) using TRIzol Universal reagent (TIANGEN, Beijing, China). The extracted RNA was then reverse transcribed with the PrimeScript RT kit (TaKaRa, Dalian, China). qRT-PCR was performed in triplicate using FastStart Universal SYBR Green Master (ROX) (Roche), with GAPDH serving as an internal control. The relative mRNA levels were calculated using the 2(-Delta Delta CT) method. The Mann-Whitney U test was employed for comparing the two groups of data due to non-normal distribution. The primers utilized for qRT-PCR are listed in Table 1.

Sections of human sinonasal mucosal tissues, embedded in paraffin at a thickness of 4 μm, were prepared and analyzed to evaluate the expression levels of the target proteins CDH3, PTHLH, PDCD4, and AR.

For immunohistochemical staining, the paraffin sections were treated with specific primary antibodies. Each slide was randomly examined and photographed in 5 high-power fields of view. Subsequently, the images were exported, and the average optical density value within the field of view was analyzed using Image-ProPlus 6.0 software, with the results being recorded. The characteristic gene expression in nasal polyps was semi-quantitatively assessed using the H-Score formula: ∑ pi × i = percentage of weak intensity cells × 1 + percentage  of  moderate  intensity   cells ×   2 + percentage  of  strong intensity  cells × 3 . In this formula, “i” indicates the classification of positive cells, where negative cells without coloring are scored as 0 points, weakly positive cells as 1 point (light yellow), moderately positive cells as 2 points (brown-yellow), and strongly positive cells as 3 points (brown-brown). “pi” represents the percentage of positive cells for each grade.

Statistical analyses were conducted using R version 4.1.1 and GraphPad Prism 8. Univariate and multivariate logistic regression analyses were employed to evaluate the diagnostic efficacy of the predictive model. All statistical tests were two-tailed, with statistical significance set at P < 0.05.

To explore the role of ARGs in CRSwNP, we conducted a systematic assessment of their expression in CRSwNP patients using the dataset GSE136825. Our analysis identified variations in the expression levels of 35 ARGs. Notably, 23 ARGs (CEACAM5, LAMB3, CDH3, MMP2, CXCR4, SERPINA1, PTHLH, MMP11, HMOX1, CCR7, HAVCR2, BCL2L15, PRKCQ, NOX4, BAX, MMP13, MMP9, TGFB1, SPIB, UBE2C, NTRK1, SEMA7A, CEACAM3) were upregulated, while 12 ARGs (LTF, PDCD4, ITGA6, CEACAM1, CPEB2, EGF, CRYAB, IL6, AR, BMP6, VTN, MIR145) showed downregulation (Figures 2A,B). Moreover, we plotted the chromosomal locations of these 35 ARGs (Figure 2C) and conducted correlation analysis to explore their interactions, revealing significant associations among specific anoikis regulatory genes (Figures 2D,E).

The CIBERSORT algorithm, based on gene expression analysis, was employed to evaluate the variation in the proportions of 22 infiltrating immune cell types in each sample. The results indicated a significant upregulation of regulatory T cells (Tregs), M2 macrophages, resting dendritic cells, resting mast cells, activated mast cells, and neutrophils in patients with CRSwNP, while plasma cells showed a downregulation (Figures 3A,B). These findings suggest a substantial impact of CRSwNP on the immune system. Additionally, correlation analysis revealed strong associations among activated dendritic cells, resting dendritic cells, M2 macrophages, activated mast cells, neutrophils, plasma cells, Tregs, and anoikis regulators (Figure 3C), indicating a significant role of ARGs in the immune infiltration alterations observed in CRSwNP patients.

To provide a more comprehensive depiction of the expression profile of ARGs in CRSwNP, we applied the consensus clustering algorithm to classify 42 CRSwNP samples based on the expression patterns of 35 ARGs. By varying the value of k from 2 to 7, we found that at k = 2, the consensus index of the cumulative distribution function (CDF) curve displayed minimal fluctuations within a narrow range, accompanied by a notably high consensus score, indicating the optimal value of k to be 2 (Figures 4A–D). Additionally, principal component analysis (PCA) revealed significant distinctions between the two clusters. As a result, we stratified the 42 CRSwNP patients into two distinct clusters.

We examined the expression of 35 DEGs between two clusters and identified distinct patterns. Cluster 1 demonstrated high expression levels of LTF, CEACAM1, CPEB2, EGF, AR, BMP6, and MIR145, whereas cluster 2 exhibited elevated levels of CEACAM5, LAMB2, PTHLH, MMP11, and BAX (see Figure 5A). Furthermore, to assess variations in immune microenvironment characteristics between the distinct anoikis-associated clusters, we analyzed differences in infiltrating immune cells and their functions. Our findings indicate that cluster 1 was enriched with resting memory CD4⁺ T cells and M1-polarized macrophages, while cluster 2 showed a higher proportion of activated dendritic cells (Figure 5B), suggesting a unique immune profile between these clusters.

Subsequently, we conducted GSVA based on GO and KEGG gene sets. The GO analysis revealed significant upregulation of pathways including the ragulator complex, early endosome lumen, glyceraldehyde-3-phosphate metabolic process, ubiquinol cytochrome C reductase activity, proteasome core complex beta subunit complex, respiratory chain complex IV, mitochondrial proton-transporting ATP synthase complex coupling factor F O, and ketosteroid monooxygenase activity in cluster 2. Conversely, pathways related to the regulation of macrophage proliferation, mesenchyme migration, saliva secretion, taurine metabolic process, epithelial cell proliferation involved in salivary, and regulation of apoptotic process in bone marrow showed downregulation in cluster 2 (Figure 5C). The KEGG analysis revealed that pathways such as oxidative phosphorylation, Parkinson’s disease, Huntington’s disease, Alzheimer’s disease, ribosome, proteasome, cardiac muscle contraction, metabolism of xenobiotics by cytochrome P450, RNA polymerase, and pentose phosphate pathway were upregulated in cluster 2. On the other hand, pathways including taurine and hypotaurine metabolism, glycine, serine and threonine metabolism, tryptophan metabolism, vascular smooth muscle contraction, beta-alanine metabolism, inositol phosphate metabolism, glycerophospholipid metabolism, homologous recombination, systemic lupus erythematosus, metabolism of xenobiotics by cytochrome P450, and drug metabolism cytochrome P450 were downregulated in cluster 2 (Figure 5D).

We utilized three algorithms to identify candidate genes associated with anoikis from a pool of 35 DEGs related to anoikis for predicting the occurrence of CRSwNP. The lasso model results indicated that 14 genes were linked to the occurrence of CRSwNP, including LTF, PDCD4, CDH3, CXCR4, SERPINA1, PTHLH, CPEB2, CRYAB, BCL2L15, AR, PRKCQ, TGFB1, SPIB, SEMA7A, and VTN (Figures 6A,B). In parallel, support vector machine (SVM) feature vectors were generated, identifying 7 genes as anoikis variables: PTHLH, PDCD4, SEMA7A, CDH3, SERPINA1, AR, and VTN (Figure 6C). Using the random forest algorithm, we identified 13 signature genes with importance scores exceeding two, including HMOX1, BCL2L15, PTHLH, CRYAB, EGF, AR, MMP9, PDCD4, NTRK1, UBE2C, BAX, BMP6, and CDH3 (Figure 6D). Finally, by intersecting the genes obtained from the three machine learning models, we identified 4 hub genes (PDCD4, CDH3, PTHLH, and AR) for further analysis (Figure 6E). The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) results indicated that these genes were predominantly involved in nuclear division, spindle organization, tubulin binding, and cell cycle regulation (Figure 6F).

Following this, we conducted receiver operating characteristic (ROC) analysis and calculated the area under the curve (AUC) values to assess the accuracy of each diagnostic gene. Our results showed that all 4 key genes demonstrated relatively high predictive values in the training set (GSE136825; Figures 7A–D). Additionally, we validated these findings in another dataset (GSE179625; Figure 7E).

A nomogram was developed for CRSwNP based on four key genes (Figure 8A). Each key gene in the nomogram is assigned a score, and the total score is calculated by summing the scores of all the characteristic genes. This total score is indicative of varying risks associated with CRSwNP. Calibration curves confirmed the accuracy of the nomogram in predicting the onset of CRSwNP (Figure 8B). Decision curve analysis illustrated the potential benefits of using the nomogram for patients with CRSwNP (Figure 8C).

The correlation between gene expression and immune cell infiltration of the four hub genes was analyzed in the combined database of GSE136825 and GSE179625. The results indicated a positive correlation of the AR gene with Macrophages M1 and a negative correlation with Macrophages M2 (Figure 9A). The CDH3 gene exhibited positive correlations with T cells regulatory (Tregs), B cells naive, Macrophages M0, and NK cells activated, and negative correlations with T cells CD8⁺, Mast cells activated, Eosinophils, and Monocytes (Figure 9B). PDCD4 gene showed a positive correlation with B cells naive, Mast cells activated, and NK cells activated, while PDCD4 gene was negatively correlated with T cells CD4⁺ memory activated, Neutrophils, and B cells memory (Figure 9C). The PTHLH gene displayed positive correlations with Macrophages M2, T cells regulatory (Tregs), T cells gamma delta, and T cells CD4⁺ memory resting, whereas the PTHLH gene showed a negative correlation with Mast cells activated (Figure 9D). These findings suggest that the expression of these genes may influence the levels of immune cell infiltration, indicating their potential involvement in immune regulation in the pathogenesis of CRSwNP.

Subsequently, we constructed the gene-miRNA regulatory network (Figure 10A) and gene-transcription factor (TF) regulatory network (Figure 10B) for the four key genes, respectively. The results revealed the intricate involvement of multiple miRNAs and TFs in regulating these pivotal genes. Furthermore, a drug interaction regulatory network for AR was developed (Figure 10C). Collectively, these findings provide valuable insights into the prospective utilization of these genes in disease diagnosis, subtype classification, survival prognosis, drug responsiveness assessment, and other pertinent domains.

Finally, the potential targets mentioned previously were validated using qRT-PCR (Figure 11A), IF (Figure 11B) and IHC (Figure 11C). As depicted in the figures, mRNA and protein levels of PTHLH and CDH3 were notably elevated in nasal polyp tissue compared to the control group, whereas PDCD4 and AR were downregulated (P < 0.05). These experimental findings align with our collected data.