Datasets GSE136825 and GSE194282, associated with CRSwNP, were retrieved from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/gds). The GSE136825 dataset, which was accessed on 17 April 2024, based on the GPL20301 platform, included inferior turbinate tissue samples from 28 healthy controls and 75 CRSwNP patient samples, 42 of which were nasal polyp tissue samples and 33 were inferior turbinate tissue samples. This dataset served as the training set for the study. Ethical approval for this study was granted by the Ethics Committee of the Affiliated Hospital of Yunnan University (Second People’s Hospital of Yunnan Province, Yunnan Eye Hospital) (Approval number: 2024229), in compliance with the Declaration of Helsinki, and oral informed consent was obtained. The GSE194282 dataset, based on the GPL17692 platform, consisted of nasal polyp tissues from 7 CRSwNP patients and hook tissue samples from 7 healthy controls (accessed on 25 April 2024), and was utilized as the validation set (Zhou et al., 2023). Additionally, 35 MP-RGs and 2,752 metabolism-related genes (MRGs) were compiled from the literature (Zhao et al., 2021; Possemato et al., 2011) and were incorporated into this study.

The DESeq2 package (v 1.38.0) (Love et al., 2014) was employed to conduct differential expression analysis on CRSwNP and normal samples in the training set, with DEGs identified based on |log₂Fold Change (FC)| > 1 and false discovery rate <0.05. Subsequently, the ggplot2 package (v 3.5.0) (Gustavsson et al., 2022) was used to generate volcano plots, and the pheatmap package (v 1.0.12) (Gu and Hübschmann, 2022) was utilized to create heat maps, highlighting the top 10 upregulated and downregulated genes sorted by log₂FC.

The scores of MP-RGs in the training set were determined using the single-sample gene set enrichment analysis algorithm from the GSVA package (v 1.46.0) (Hänzelmann et al., 2013). Significant differences in these scores between CRSwNP and normal samples (P < 0.05) prompted the use of WGCNA through the WGCNA package (v 1.71) (Langfelder and Horvath, 2008) to analyze expression values associated with the MP-RGs scores. Hierarchical clustering was performed using Euclidean distances to detect potential outliers, which were subsequently removed. The MP-RGs scores were then re-clustered as sample traits. To construct the scale-free network, the best soft threshold was determined by ensuring that the scale-free fit index (R ²) was 0.9 and that the average connectivity was nearly zero. The minimum number of genes in each gene module was set to 100, according to the rules of the hybrid dynamic tree cutting algorithm. Gene similarity dendrograms were constructed to calculate the relationship between module feature genes and traits, and heatmaps were generated to visually depict these correlations. Finally, modules with a correlation to MP-RGs scores exceeding 0.4 and a p-value of less than 0.05 were selected as key modules. The genes within these key modules were identified as the key modular genes linked to MP-RGs.

Using the ggvenn package (v 0.1.9) (Zheng et al., 2022), the intersection of DEGs, MRGs, and key modular genes was identified. These intersected genes were subsequently utilized as candidate genes for further analysis.

Following the identification of the candidate genes, clusterProfiler software (v 4.7.1.003) (Yu et al., 2012) was used to perform Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The GO analysis assessed biological processes (BP), cellular components, and molecular functions (MF). The top 10 enriched pathways were selected and displayed based on the results. To explore the interactions of the candidate genes at the protein level, the STRING database (http://string-db.org) (accessed on 17 April 2024) was utilized. This facilitated the construction of a protein-protein interaction (PPI) network, in which only interactions with a score above 0.4 were considered. The network was then visualized using Cytoscape software (v 3.8.2) (Liu et al., 2020). The Cytoscape plugin CytoNCA was employed to identify the most significant nodes in the PPI network. This process involved applying a range of parameters, including degree centrality, betweenness centrality, proximity centrality, eigenvector centrality, local average connectivity-based approach, and network centrality. The criterion established was that the threshold for each centrality metric was greater than the median. After applying these criteria, a PPI network containing the core genes was obtained.

To identify genes closely associated with CRSwNP, a least absolute shrinkage and selection operator (LASSO) regression method penalized regression model was constructed using the glmnet package (v 4.1.4) (Li et al., 2022) based on the training set. With “CRSwNP disease status” as the dependent variable and the expression levels of core genes as independent variables, the regularization parameter λ was optimized through 10-fold cross-validation. Genes unrelated to CRSwNP had their coefficients shrunk (with coefficients approaching zero), while genes making significant contributions to classification were retained, thereby reducing the complexity of the model. Ultimately, genes with stable predictive ability were identified. Similarly, with “CRSwNP disease status” as the prediction phenotype, the randomForest package (v 4.7.1.1) (Alderden et al., 2018) was employed to carry out a random forest (RF) analysis. This analysis was evaluated using out-of-bag (OOB) error. When the number of decision trees reached approximately 600, the error stabilized, and the optimal number of trees was determined accordingly. The top 10 genes ranked by importance were selected for subsequent analyses. Additionally, to prevent overfitting, multiple decision trees were constructed by bootstrap sampling of samples and random selection of features, aiming to reduce the bias of individual trees. Finally, the ggvenn package was used to identify the signature genes for this study by analyzing the crossover results.

To evaluate the ability of the signature genes to distinguish between CRSwNP and normal samples, the pROC package (v 1.74.0) (Robin et al., 2011) was employed to plot receiver operating characteristic (ROC) curves for these genes in both the training and validation sets. The area under the curve (AUC) for these ROC curves was subsequently calculated. Additionally, expression validation of signature genes was conducted in both the training and validation sets. Signature genes with an AUC >0.7, consistent and significantly different expression patterns across all datasets, were considered potential biomarkers. A nomogram was then constructed based on these biomarkers using the rms package (v 6.5.0) (Sachs, 2017). Each biomarker was assigned a specific score within the model, with the cumulative score of all biomarkers contributing to the total score. To assess the nomogram’s predictive performance, calibration curves were constructed. An ideal model would exhibit calibration curves with a slope approaching 1, indicating high accuracy in prediction.

To determine the biological pathways and functions of biomarkers in the training set, the psych package (v 2.2.9) (Author anonymous, 2023) was employed to evaluate the correlation between the biomarkers and other genes and to rank the correlation coefficients. KEGG gene annotations from the R package “org.Hs.eg.db” (v 3.16.0) (Qing et al., 2022) were used as the background, and the genes were matched to this set. Subsequently, GSEA was performed on the biomarkers using the clusterProfiler package to identify pathways that were significantly enriched (P < 0.05).

To explore the molecular regulatory mechanisms associated with the biomarkers, the miRDB (https://www.mirdb.org) and miRWalk (http://mirwalk.umm.uni-heidelberg.de/) databases (accessed on 17 April 2024) were utilized to predict miRNAs that interact with these biomarkers. The intersection of the predicted results from both databases was considered the set of key miRNAs. Based on the predicted key miRNAs, a miRNA-mRNA network was constructed. Subsequently, StarBase (http://starbase.sysu.edu.cn/index.php) (accessed on 17 April 2024) was used to predict the long non-coding RNAs (lncRNAs) that function upstream of these key miRNAs. An lncRNA-miRNA interaction network was then built to map the regulatory relationships. The two networks were subsequently integrated to form the lncRNA-miRNA-mRNA network. Additionally, to explore transcription factors (TFs) that target regulatory biomarkers, the JASPAR database (http://www.jaspar.genereg.net) (accessed on 17 April 2024) was used to identify TFs that regulate the biomarkers, and a TF-mRNA network was constructed. The results of these analyses were graphically represented using Cytoscape software.

To obtain biomarker-disease associations, biomarkers were imported into the Disease Gene Network database (http://www.disgenet.org/) (accessed on 17 April 2024) and analyzed, with a score >10 indicating an association. Additionally, to identify potential drugs targeting the biomarkers, relevant drug compounds were predicted using the Drug Gene Interaction Database (accessed on 17 April 2024). Subsequently, a drug-biomarker interaction network was delineated using Cytoscape software to visually represent the relationships between the drugs and biomarkers.

The expression of these biomarkers was confirmed using RT-qPCR. Five pairs of CRSwNP and control tissue samples were obtained from the Affiliated Hospital of Yunnan University. Informed consent was obtained from all participants, and the study was approved by the hospital’s ethics committee (approval number: 2024229). Total RNA was extracted from 10 samples using TRIzol reagent (Ambion, Austin, United States) following the manufacturer’s guidelines. The quality and concentration of extracted RNA were verified using a NanoPhotometer N50, with OD260/280 ratios ranging from 1.8 to 2.0, indicating high RNA purity. Next, first-strand cDNA was synthesized from 1 μg of total RNA using the SweScript First-Strand cDNA Synthesis Kit (Servicebio, Wuhan, China) according to the manufacturer’s protocol: the reaction mixture (20 μL total volume) included 4 μL of 5× reaction buffer, 1 μL of reverse transcriptase, 1 μL of oligo (dT) primer, 1 μg of RNA template, and nuclease-free water to volume. RT-qPCR was performed using the 2 × Universal Blue SYBR Green qPCR Master Mix (Yeasen, China) on a CFX96 real - time PCR detection system (BIO - RAD, United States). The 10 μL reaction system contained 5 μL of 2xUniversal Blue SYBR Green qPCR Master Mix, 1 μL of forward primer (10 μM), 1 μL of reverse primer (10 μM), 3 μL of cDNA template. The thermal cycling conditions were: initial denaturation at 95 °C for 1 min, followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 55 °C for 20 s, and extension at 72 °C for 30 s. After the reaction was completed, specificity was verified via melting curve analysis. The specific primer sequences for target biomarkers and the internal reference gene are listed in Table 1. All primers were designed using NCBI Primer-BLAST. During design, the primers were ensured to span exon-intron boundaries to exclude genomic DNA contamination. Besides, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the internal reference gene for data normalization. Finally, the relative expression levels of target biomarkers were calculated using the 2^(−ΔΔCt) method (Livak and Schmittgen, 2001), where ΔCt = Ct (target gene) - Ct (GAPDH), and ΔΔCt = ΔCt (CRSwNP group) - ΔCt (control group). Results are presented as fold changes relative to the control group.

Statistical analyses for the bioinformatics component were performed using R software (version 4.2.2). For comparisons between two groups, the Wilcoxon rank-sum test was applied. For RT-qPCR data, relative expression levels were calculated using the 2^−ΔΔCt method and compared between groups using an unpaired t-testthe. In all statistical analyses, a two-tailed P-value <0.05 was considered statistically significant.

Differential expression analysis conducted on the training set identified 4,312 DEGs, of which 2,672 were upregulated and 1,640 downregulated (Figures 1A,B). Further assessment of the MP-RGs scores between the CRSwNP and normal samples revealed significantly higher MP-RGs scores in the CRSwNP samples compared to the normal controls (P < 0.0001) (Figure 2A). This finding suggests that macrophage polarization may play a critical role in the pathophysiology of CRSwNP. Consequently, WGCNA was employed to identify key modular genes associated with MP-RGs. As shown in the figure, no significant outliers were detected, and all samples were retained for subsequent analysis (Figure 2B). The optimal soft-threshold power was determined to be 9, where R ² reached 0.9, and the average connectivity value approached zero (Figure 2C). A dendrogram was then constructed, leading to the identification of 19 modules (Figure 2D). The MEblack module (cor = 0.49), MEred module (cor = −0.56), and MEbrown module (cor = 0.69), which exhibited absolute correlation coefficients greater than 0.4 with the MP-RGs scores, were designated as key modules. These modules collectively comprised a total of 1,693 key modular genes (Figure 2E). Metabolic reprogramming was a key regulatory mechanism underlying macrophage polarization. Therefore, after screening 4,312 DEGs and 1,693 key modular genes, 2,752 MRGs were introduced for overlap analysis, with the aim of further focusing on candidate genes that were also involved in metabolic processes from the previously screened genes, we identified 133 candidate genes (Figure 2F).

Among the 289 enriched GO-BP entries, significant associations were observed between candidate genes and functions such as “organic anion transport”, “small molecule catabolic process”, and “glycolipid metabolic process”. In the 28 enriched GO-CC entries, candidate genes were associated with structures such as the “ion channel complex”, “transmembrane transporter complex”, and “phosphatidylinositol 3-kinase complex”, among others. Within the 132 GO-MF categories, associations were found between the candidate genes and functions such as “metal ion transmembrane transporter activity”, “channel activity”, and “lipase activity”, among others (Figure 3A), indicating that the candidate genes might have acted by regulating the clearance of inflammatory mediators and the metabolic reprogramming of macrophages, thereby maintaining the chronic inflammatory microenvironment of the nasal mucosa; meanwhile, they might have synergistically participated in the signal transduction of immune cell activation and lipid metabolism disorders, promoting the abnormal proliferation of nasal polyp tissues and the persistence of inflammation. Additionally, these candidate genes were found to be enriched in 41 KEGG pathways, participating in various BPs, including metabolism, signaling, neurotransmission, and hormone secretion (Figure 3B), suggesting that candidate genes may act by participating in the crosstalk regulation of the nasal mucosal neuro-immune-endocrine-metabolic network, synergistically promoting the abnormal proliferation of nasal polyp tissues and the chronic persistence of inflammation. Furthermore, the PPI network analysis revealed interactions between the 133 candidate genes (Figure 3C). Subsequently, a total of 55 core genes were identified through screening with the CytoNCA plugin in Cytoscape (Figure 3C).

Based on the 55 core genes, signature genes were further screened using LASSO regression and RF machine learning algorithms, respectively. A total of 11 genes, including PLA2G2A, NOS3, KCNAB1, ALOX5, PIK3C2G, HMOX1, CA4, PLA2G7, GNMT, BCAT1, and ITPKB, were selected from the 55 core genes by the LASSO regression algorithm, with lambda. min (0.09792358) as the optimal model (Figures 4A,B). The decision tree exhibiting the lowest error rate was derived using RF, and the importance of the genes was ranked. The top 10 genes identified were HMOX1, PLA2G7, SLC12A2, ADCY7, GNMT, PLA2G2A, ALOX5, KCNJ8, ALOX5AP, and MGLL (Figures 4C,D). Finally, the results from both algorithms were intersected to identify five signature genes: PLA2G2A, ALOX5, HMOX1, PLA2G7, and GNMT (Figure 4E).

The diagnostic value of the five signature genes was further analyzed and evaluated. The results demonstrated that ALOX5, HMOX1, and PLA2G7 exhibited AUC values of 0.818, 0.902, and 0.853, respectively. These values exceeded 0.7, indicating that these three genes possessed a strong capacity to differentiate between CRSwNP and normal samples. This is depicted in Figure 5A. Similarly, the AUCs of these three signature genes in the validation set surpassed 0.7, further corroborating their robust diagnostic potential (Figure 5B). To refine the identification of biomarkers, the expression levels of the five signature genes were assessed in both the training and validation sets. Notably, the expression of these three genes was significantly higher in the CRSwNP group compared to the normal group in both datasets (P < 0.05) (Figures 5C,D). This observation highlights the promise of these genes as diagnostic markers for CRSwNP, given their consistent and pronounced expression differences between the CRSwNP and normal groups across both datasets. Consequently, ALOX5, HMOX1, and PLA2G7 were identified as the biomarkers in this study.

Subsequently, a nomogram for CRSwNP was constructed based on ALOX5, HMOX1, and PLA2G7 (Figure 6A). The calibration curve of the nomogram exhibited a slope close to 1, demonstrating that the prediction accuracy of the nomogram was excellent (Figure 6B).

GSEA analysis revealed that the biomarkers ALOX5, HMOX1, and PLA2G7 were significantly associated with 1,859, 2,228, and 2,025 pathways, respectively, in the training set. The top 10 pathways for each biomarker are presented in Figures 7A–C. Upon examining the significantly enriched pathways for all biomarkers (Supplementary Tables S1-S3), it was observed that they co-enriched into immune cell-related pathways, including ‘immunoregulatory interactions between a lymphoid and a non-lymphoid cell’ and ‘natural killer cell-mediated cytotoxicity’, among others. It was concluded that the three biomarkers may be strongly associated with inflammation and influence immunological pathways that impact the development and progression of CRSwNP.

Two databases were queried to predict key miRNAs potentially targeting the biomarkers, resulting in the identification of 10 key miRNAs corresponding to the three biomarkers. The miRNA-mRNA regulatory network revealed that no miRNAs targeted all biomarkers simultaneously (Supplementary Figure S1A). Subsequent predictions of the upstream lncRNAs for the 10 miRNAs identified that only hsa-miR-642a-5p predicted 10 lncRNAs, including FOXC1 and NORAD (Supplementary Figure S1B). A lncRNA-miRNA-mRNA regulatory network containing 22 nodes and 20 edges was further constructed, visualizing complex interactions such as ZFAS1-hsa-miR-642a-5p-PLA2G7 (Figure 8A). Additionally, the JASPAR database was utilized to identify TFs for the regulatory biomarkers, leading to the construction of a TF-mRNA regulatory network (Figure 8B). Among these, ALOX5 and HMOX1 were associated with the most TFs, and FOXC1 simultaneously regulated all three biomarkers.

Analysis of the biomarker-disease network map revealed that ALOX5, HMOX1, and PLA2G7 were strongly associated with asthma (Figure 8C). Furthermore, a drug-biomarker interaction network was constructed by screening drugs from the DrugBank database. Examination of this network demonstrated that ‘selenium’ was identified as a co-target of both ALOX5 and HMOX1 (Figure 8D).

The expression of biomarkers in clinical samples was further analyzed using RT-qPCR. Melting curve analysis revealed a single, sharp peak for ALOX5, HMOX1, PLA2G7 and the internal reference gene GAPDH, indicating specific amplification without non-specific products or primer-dimers (Supplementary Figure S2). The results demonstrated that the expression levels of ALOX5, HMOX1, and PLA2G7 were significantly higher in CRSwNP samples compared to control samples (P < 0.05) (Figure 9). These findings from the clinical sample validation are consistent with the results of the bioinformatics analysis, thereby enhancing the reliability of the bioinformatics approach.