Introduction

Chronic rhinosinusitis with nasal polyposis (CRSwNP) represents a prevalent localized persistent inflammatory disorder. Patients afflicted with this condition commonly experience symptoms featuring nasal blockage, rhinorrhea, olfactory dysfunction, and pain in the face, which markedly impair their quality of life and work efficiency. Additionally, the disease incurs substantial economic burdens.^1,2 Over the past years, the occurrence of CRSwNP has escalated significantly due to alterations in people’s living environments and lifestyles, posing a grave threat to their physical and mental well-being.³ The estimated global prevalence of chronic rhinosinusitis(CRS) spans from 5% to 12%.⁴ Genetic analysis of CRSwNP can facilitate the elucidation of genes implicated in modulating the disease process, thereby enabling the development of more precise therapeutic approaches and biomarkers.⁵ Considering the high prevalence of CRSwNP, the severe discomfort experienced by patients, and the pressing need for treatment, it is urgent to detect efficacious diagnostic predictors and therapeutical targets.

High-throughput sequencing advances have led to increased research combining sequencing and bioinformatics. These tools help identify key genes and pathways in diseases or biological processes, providing a basis for early diagnosis and drug development.In previous bioinformatics analyses, genes such as XIST, TAS2R19, TYROBP, and MAP1B⁶-⁹ have been identified to facilitate the progression of CRSwNP. Simultaneously, the application of supervised or unsupervised machine learning algorithms has demonstrated their potential for unveiling concealed relationships within high-dimensional data.^10,11 Moreover, machine learning holds significant value in the evaluation of high-dimensional transcriptome data and the ascertainment of the positions of genes with biological relevance.^12,13 In the field of CRS, it has been applied to research on the burden of medical care¹⁴ and the prediction of olfactory loss.¹⁵ Nonetheless, it has yet to be utilized in the mechanistic investigations of CRSwNP.

This study is predicated on transcriptome data derived from the Gene Expression Omnibus (GEO). Bioinformatics methodologies were employed to detect key genes implicated in CRSwNP, and Weighted gene co-expression network analysis (WGCNA) was utilized to pinpoint crucial gene modules linked to CRSwNP phenotypes. The outcomes of the screened differentially expressed genes (DEGs) and WGCNA underwent gene ontology (GO) enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, and disease ontology (DO) enrichment analysis. Different results will be obtained when analyzed with different logical algorithms. In order to compensate for the impact of the shortcomings of different algorithms on the analysis results as much as possible, three representative machine learning algorithms, namely, least absolute shrinkage and selection operator (LASSO), randomForest, and support vector machine recursive feature elimination (SVM-RFE), were employed to ascertain feature genes, which were subsequently validated for diagnostic efficacy through ROC analysis and GSE136825 data.

Murine models, which provides a bridge between genetic predictions and experimental validation, are widely used in the research of CRSwNP,^16–18so a murine model was constructed, and qPCR experiments were conducted to validate the feature genes, which demonstrated that BTBD10, ERAP1, GIPC1, and PEX6 may serve as key genes in CRSwNP. The ssGSEA algorithm was leveraged for analyzing the expression level exhibited by immune cell genes, successfully identifying and quantifying the diverse immune cells residing in the nasal mucosa of CRSwNP and control samples, and characterizing the link between immune infiltration and CRSwNP biomarkers.

While previous studies have indeed confirmed a number of potential target genes related to CRSwNP, the precise genes and that are intricately associated with the condition are yet to be fully elucidated. The intention of this study is to detect prospective biomarkers for CRSwNP by screening key genes and signaling pathways pertinent to the disease employing bioinformatics techniques. Furthermore, we have carried out validation. This will facilitate the advancement of CRSwNP treatment and mechanistic research. The study’s flowchart is illustrated in Figure 1.

Figure 1 Flowchart of the study process GO, Gene Ontology; DO, Disease Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; SVM-RFE, support vector machine recursive feature elimination; RF, random forest; ROC, receiver operating characteristic.

Materials and Methods

Data Download

CRSwNP datasets GSE136825 and GSE107624 both were procured from the NCBI’s GEO (https://www.ncbi.nlm.nih.gov/geo/). The GSE136825 dataset encompasses 42 nasal polyposis (NP) samples and 28 normal control samples, whereas the GSE107624 dataset contains 21 NP samples and 21 normal control samples. In the present study, GSE107624 was utilized as a training set, while GSE136825 functioned as a validation set. All datasets underwent standardized data preprocessing.

Differential Expression Analysis

The differentially expressed genes (DEGs) within GSE107624 were identified through the employment of the “Limma” R package, with the screening criteria established as |log2 fold change (FC)| > 0.5 and p < 0.05. Heatmaps and volcano plots of the DEGs were created by employing the “Pheatmap” R package and the “ggplot2” R package, respectively.

Weighted Correlation Analysis

The “WGCNA” software package in R was utilized to examine gene expression trends among various samples and to identify potential biomarkers or therapeutic targets. Initially, the existence of each gene was investigated, and an adjacency matrix was devised to measure the degree of correlation strength between nodes. Thereafter, the adjacency matrix was converted into a topological overlap matrix (TOM) to numerically the likeness within nodes. The dissimilarity scores were computed to facilitate hierarchical clustering analysis. Utilizing a dynamic tree cutting approach with a threshold module size of 50, we identified co-expressed gene modules. Subsequently, we assessed the association between these gene modules and CRSwNP by evaluating gene significance (GS) and module membership (MM) values. In the end, the module demonstrating the strongest correlation with the disease was designated as the key module.

Enrichment Analysis of DEGs

The DOSE package in R was adopted for exploring biological functions in GO, DO, and KEGG enrichment analyses. Gene Set Enrichment Analysis (GSEA) was employed to delve into the modifications in biological functions and pathways between CRSwNP and normal control samples. The limit for statistical significance was set at a P-value < 0.05.

Selection of Core Feature Genes (Hub Genes)

Three machine learning algorithms, specifically LASSO, SVM-RFE, and randomForest, were utilized to screen for hub genes in CRSwNP. The “glmnet” package was utilized to apply the LASSO algorithm, known for its efficient execution of LASSO and elastic-net regularization techniques, set the α value of the glmnet function to 1, obtained the best λ by ten cross-validations finally obtained the hub genes based on the best λ. For the SVM algorithm, the “e1071” package was employed, offering a wide array of utilities for SVM model training and classification operations. The RF algorithm was executed with the help of the “randomForest” package. An error rate assessment was conducted for each tree, with the optimal trees being identified as the count that resulted in the minimum error rate while ensuring model robustness. Subsequently, the genes screened out by these machine learning algorithms were overlapped to identify potential hub genes for further investigation. The expression rates of the hub genes were scrutinized in GSE107624.

Assessment of the Diagnostic Efficacy of Hub Genes

A nomogram was crafted using the expression levels of hub genes, and the risk of developing CRSwNP was predicted via the “rms” package in the R software. A calibration plot was constructed to assess the predictive precision of the model. Subsequently, curves for decision analysis and clinical impact were drawn to determine whether the application of the diagnostic model would result in beneficial clinical decisions for patients. The diagnostic merit of hub genes and the nomogram was then evaluated by establishing ROC curves, with the area under the curve (AUC) and 95% CI worked out.

Examination of Immune Infiltration of Hub Genes

The Single-sample gene set enrichment analysis (ssGSEA), which is an expansion of the traditional GSEA approach, has become a prevalent technique in bioinformatics research for investigating immune infiltration, implemented in the R language, was adopted to quantify the proportions of 28 types of immune cell within the merged dataset. Violin plots were employed to display the contrast in immune cell composition between CRSwNP and healthy samples. The relationship between hub genes and infiltrating immune cells was scrutinized via the R software. The findings of the analysis were presented in a visual format via the “ggpubr” package in R.

Animals

Male C57BL/6J mice, aged 6 to 8 weeks and free from specific pathogens, were supplied by Jinan Pengyue Experimental Animal Breeding Co., Ltd. and maintained at the Experimental Animal Center of Yuhuangding Hospital. The mice were allowed to adapt to the animal facility environment for a period of one week prior to the commencement of the experiment, which had been authorized by the Ethics Committee of Yantai Yuhuangding Hospital, affiliated with Qingdao University. All animal experiments were conducted in accordance with the ARRIVE guidelines, the UK Animals (Scientific Procedures) Act 1986 and related guidelines, and the EU Directive on Animal Experiments 2010/63/EU. The mice were housed in temperature-controlled environments ranging from 22°C to 24°C, with a 12-hour light/dark cycle from 7:00 am to 7:00 pm, and they had free access to food and water.

Establishment of CRSwNP Model in Mice

With regard to the meticulous experimental setup, the techniques utilized in prior research were replicated.¹⁶ Six-week-old mice in the NP group were intranasally challenged with a mixture of Aspergillus proteinase (AP) and ovalbumin (OVA). A solution containing 2 units of AP (Sigma-Aldrich, St. Louis, MO) and 75 mg of OVA (Sigma-Aldrich) was prepared by diluting it in sterile PBS to a total of 20 mL and was then administered intranasally to each mouse on three occasions per week for a duration of 12 weeks, with 6 mice in each experimental group. The mice in the normal group (6 mice per group) were dosed with PBS at equal intervals. Following nasal instillation, the mice were held with their heads downward to prevent the reagents from entering the lungs.

Sample Collection

Mice were induced into anesthesia with a 2% dose of isoflurane and subsequently euthanized. Nasal mucosa samples were collected, and the dermal tissue from the mouse’s head was separated by blunt incision under a microscope to expose the mouse’s skull. The connection between the zygomatic bone and the maxillary structures on either side of the skull was then trimmed with scissors. The nasal bones above the nostrils were separated with forceps, and the soft tissue covering the palate was also separated with forceps to expose the bottom of the nasal cavity. Additionally, the alveolar bone of the maxilla was severed with forceps. The upper jawbone and the nasal tip were removed with scissors. At the base of the maxilla, a blade or a pair of scissors was employed to cut the bone septum free from the maxillae on the left and right. Forceps were utilized to grasp the front end of the maxilla and pry the maxilla on both sides. They were also employed to grip the base of the nasal septum and detach the nasal septum. It is of paramount importance to exercise caution to avoid damaging the nasal mucosa. The mucosa lining the nasal cavity over the maxillae on either side was observable. Fine forceps were employed to collect and preserve the nasal mucosa on both sides of the bony nasal septum. The complete array of tissue samples was stored in a −80°C freezer for further research.

Histological Analysis

A 4% PFA solution was used to fix the nasal mucosa tissue of mice for a duration of 12 hours, and then a 5% nitric acid solution was employed for decalcifying the material for a duration of three days. The tissue was embedded in paraffin in line with standard practices. A paraffin slicer (RM2016, Leica, Germany) was used to cut 4-μm thick sections. The sections underwent hematoxylin and eosin (H&E) staining to facilitate the evaluation of polypoid lesions. Histological analysis of H&E was performed by two impartial pathologists unaware of the study design.

Total RNA Extraction from Mouse and Quantitative Real-Time Analysis of PCR Expression

The RNeasy mini-kit (SparkJade, China) was employed to isolate total RNA from mouse olfactory bulbs (OB). cDNA was subsequently produced via the QuantiTect Reverse Transcription Kit (SparkJade, China). Real-time PCR quantification was carried out with SYBR Green Master Mix (SparkJade, China) on the QuantStudio 3 Real-Time PCR System. (SparkJade, China). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) served as the internal control benchmark. Listed below were specific primers (Sangon Biotech, Shanghai, China):

Mouse PEX6: F: CACACTCCGTCATCTCCTCCTTG; R: GACAGCACTGGCACCTCTAGC

Mouse BTBD10: F: AGGAACAGCAGCCAGTCAAGC; R: GTCACTCGCTCAGATGTCCTCAC

Mouse CYFIP1: F: ACACGCTACAACTACACCACTGAG; R: CACACTCTCCATCCTGCCCATC

Mouse ERAP1: F: GCCAACACTCATCATCAACCTCAC; R: TGCACATTCCGCCCACTCAC

Mouse GIPC1: F: CCTCGCCTCGTGTTCCATACC; R: ATCTTGCCGTACAGCTCCTTGAC

Mouse NYNRIN: F: TCTGGCTATTGCTTCTACCGTGAC; R: AGGTGGGCGTATGTGGTTGTG

Analysis of Statistical Data

Comprehensive statistical analyses and graphical illustrations were undertaken via R software (version 4.1.2). Data analysis was conducted using GraphPad Prism 8.0 (GraphPad Software, La Jolla, CA). The diagnostic performance of hub genes was assessed by employing ROC curve analysis. Unpaired t-tests were conducted to evaluate the differential expression of core genes. Unless otherwise stated, P < 0.05 was defined as statistically significant.

Results

Selection of DEGs and GSEA

A total of 239 DEGs were obtained from GSE107624, with 170 upregulated and 69 downregulated in CRSwNP. Volcano plots and heatmaps were constructed (Figure 2A and B). GSEA demonstrated that DEGs exhibiting upregulation might be predominantly in association with aminoacyl-tRNA biosynthesis, antifolate resistance, nucleocytoplasmic transport, pyrimidine metabolism, and ribosome biogenesis in eukaryotes. The downregulated DEGs might be chiefly correlated with the metabolic processes of drugs − cytochrome P450, the metabolism of xenobiotics by cytochrome P450, nitrogen metabolism, and pentose phosphate pathway (Figure 3A). GSEA ridge plots were also constructed (Figure 3B).

Figure 2 The volcano plot and the heatmap of DEGs screened from GSE107624. (A) The heatmap representing the differentially expressed genes (DEGs). The red coloration denotes genes with high expression, and the blue coloration denotes genes with low expression. (B) The volcano plot is a graphical illustration of the differentially expressed genes. The red coloration represents genes with high expression, and the green coloration represents genes with low expression.

Figure 3 GSEA as a Tool for Analyzing DEGs. (A) The GSEA plot presenting the enriched signaling pathways based on the hallmark gene sets. (B) The ridge plot of the top 30 most enriched pathways in CRSwNP patients.

WGCNA Analysis and Detection of Modules

All the genes were ordered in a descending manner of variance, and the leading 25% of genes in terms of variance were selected for analysis. Following the removal of outliers, a sample clustering tree (Figure 4A) was generated. A soft threshold of 7 (R² = 0.85) was set to fabricate a network with a scale-free nature. Subsequently, an adjacency matrix was constructed and a matrix reflecting topological overlap (Figure 4B) was generated. Functionally similar modules were determined through the combination of average hierarchical clustering and dynamic tree clipping (Figure 4C). The magenta module exhibited the highest correlation of module eigengene (ME) with CRSwNP (r =0.44; p =0.004), followed by the darkturquoise module (r = 0.4; p = 0.008, Figure 4D). Consequently, the status of these two modules as key modules for CRSwNP has been acknowledged. Furthermore, the linkage between module membership (MM) and gene significance (GS) was calculated within the magenta and darkturquoise modules (Figure 4E and F).

Figure 4 Establishment of WGCNA and Detection of Key Modules. (A) Dendrogram of sample clustering. (B) Examination of the scale-free index and average connectivity across a range of soft threshold powers (β). (C) Dendrogram representing the grouping of all differentially expressed genes, characterized by measurement variability (1-TOM). The colored ribbons represent the findings derived from the automated single-block analysis. (D) A heatmap illustrating the correlation between modules and CRSwNP. The numbers within and outside the parentheses demonstrate the p-values and correlation coefficients in respective order. The linkage between MM (X-axis) and GS (Y-axis) for genes within the magenta module (E) and the darkturquoise module (F).

DO, GO, and KEGG Enrichment Analysis of Overlapped Genes

The overlap of DEGs with genes from the magenta and darkturquoise modules, identified through WGCNA, yielded 56 common genes (Figure 5A). The 56 overlapped genes were subject to DO, GO, and KEGG functional enrichment analysis via the “clusterProfiler” package in R. The DO enrichment analysis revealed that the overlapped genes might be related to conditions including fatty liver disease, Parkinsonism and nephroblastoma (Figure 5B). The GO enrichment analysis established that these genes were predominantly enriched in such biological processes (BP) as nucleus organization, and the regulation of DNA methylation. In terms of cellular components (CC), these genes were chiefly enriched in cytoplasmic ribonucleoprotein granule, and nuclear periphery, among others. Regarding the molecular functions (MF), these genes were characterized by the enrichment in amide binding, and peptide binding (Figure 5C). KEGG enrichment analysis revealed their associations with pathways such as Nucleocytoplasmic transport, Dopaminergic synapse, and HIF-1 signaling pathway (Figure 5D).

Figure 5 DO, GO, and KEGG Enrichment Analysis of Overlapped Genes. (A) Venn diagram illustrates the overlap of DEGs with Turquoise module-specific genes. (B) The findings of the DO enrichment analysis. (C) The outcomes of the GO enrichment analysis. (D) The outputs of the KEGG enrichment analysis.

Screening of Hub Genes Based on Three Different Machine Learning Algorithms

Initially, SVM-RFE was employed to screen the aforementioned 56 differential feature genes. The minimum 5-fold cross-validation error resulted in the detection of 21 feature genes (Figure 6A). Secondly, the random forest algorithm identified 40 feature genes (Figure 6B). Finally, LASSO regression analysis indicated that 13 feature genes were included in the regression model when selecting the log λ with the minimum mean square error (Figure 6C). The genes identified by SVM-RFE, random forest, and LASSO analyses exhibited significant overlap. Finally, a total of 11 hub genes were identified in GSE107624 (Figure 6D), namely BTBD10, RAB6A, SPEN, TMEM87B, COBLL1, CYFIP1, ERAP1, GIPC1, NYNRIN, PEX6 and VHL (Figure 6E). All of these were of statistical significance.

Figure 6 Screening of Hub Genes Using Machine Learning Algorithms. (A) SVM-RFE analysis. (B) Random forest analysis. (C) LASSO regression analysis. (D) A Venn diagram illustrating the overlapping analysis through SVM-RFE, random forest, and LASSO. (E) A bar chart presenting the differential expression of hub genes in GSE107624. NP: nasal polyposis.

Construction of a CRSwNP Prediction Model and Appraisal of the Diagnostic Reliability of Key Genes and the Prediction Model

A nomogram was built upon the 11 hub genes using logistic regression analysis (Figure 7A). The calibration plot demonstrated a marked degree of concordance between the calibration prediction curve and the standard curve (Figure 7B). The decision curve analysis (DCA) indicated that the model had a favorable clinical benefit (Figure 7C). In GSE107624, the AUCs values of the prediction model for 11 hub genes were also >0.7, indicating a desirable predictive performance (Figure 7D). In addition, the diagnostic efficacy of these key genes and the prediction model in the validation set GSE136825 was evaluated (Supplementary Figure 1). These phenomena indicated that the screened signature genes had remarkable diagnostic efficiency in forecasting CRSwNP.

Figure 7 Construction and Validation of the Prediction Model. (A) The nomogram model built upon hub genes in GSE107624. (B) Calibration curve of the nomogram. (C) DCA. (D) Evaluation of the diagnostic efficacy of key genes and the prediction model in GSE107624. DCA: decision curve analysis.

Results of Infiltrative Behavior of Immune Cell

The ssGSEA was employed to measure the extent of immune infiltration. In comparison to the control group, the CRSwNP group showed decreased infiltration levels of Type 2 T helper cell, Activated CD4 T cell, Memory B cell, Mast cell, and Natural killer T cell. Furthermore, the immature dendritic cell exhibited a association with the infiltration levels of Effector memory CD8 T cell (Figure 8A). Examination of correlation revealed that there were affirmative correlations between CYFIP1 and type 2 T helper cell. CYFIP1 and BTBD10 were positively correlated with memory B cell. ERAP1 presented a positive connection with natural killer T cell; TMEM87B demonstrated a negative correlation with Immature dendritic cell. NYNRIN, COBLL1, GIPC1, SPEN, CYFIP1, and BTBD10 exhibited a negative association with Effector memory CD8 T cell (Figure 8B).

Figure 8 Results of Immune Infiltration. (A) Infiltration rates of 28 immune cell categories in the non-diseased group versus that in the CRSwNP group. (B) A heatmap illustrating the linkage of target genes to immune cells. *: p<0.05; **: p<0.01; ***: p<0.001; ****: p<0.0001; ns: non-significant.

Expression Levels of Hub Genes in CRSwNP Mouse Model

To validate the observed changes in gene expression GSE107624, we utilized GSE136825 as a validation set. It was observed that in both datasets, the gene expression of BTBD10, CYFIP1, ERAP1, GIPC1, NYNRIN and PEX6 exhibited an similar upward trends (Supplementary Figure 2). To verify the expression profiles of hub genes in CRSwNP, a mouse experimental model was established. Examination of CRSwNP mouse nasal tissues via H&E staining revealed the presence of polypoid lesions and diffuse mucosal swelling, whereas normal mice exhibited a smooth thin mucosa (Figure 9A). Nasal mucosal samples were obtained from nine control mice and nine mice with CRSwNP. The RT-qPCR findings indicated that the mRNA expression levels of BTBD10, GIPC1, NYNRIN, and PEX6 in the CRSwNP group were relatively higher than those observed in the control group (Figure 9B), which was consistent with previous findings. No significant trends were observed for CYFIP1 and ERAP1.

Figure 9 (A) H&E staining to assess the nasal mucosal status in different groups. The yellow line frames indicate the same magnified areas. The black arrows indicate the presence of uneven, swollen, and ruptured nasal mucosa. (B) The magnitude of expression for BTBD10, GIPC1, NYNRIN, and PEX6 in the healthy group and the CRSwNP group of mice. *: p<0.05; ***: p<0.001.

Discussion

NP is defined as a benign protrusion of the inflammatory nasal and sinus mucosa in CRSwNP. CRSwNP is a condition marked by the persistent inflammatory process and the emergence of edematous polyps in the nasal mucosa. The symptoms associated with this condition include nasal congestion, olfactory dysfunction, and rhinorrhea. CRSwNP influences a substantial fraction of the global community, with its incidence rising annually. This imposes a substantial burden on human health and the social economy, and has become an increasingly severe global public health issue.^1,2,19,20 Consequently, the detection of hub genes and pathways associated with the disease, as well as the exploration of new diagnostic or therapeutic targets, are of paramount importance for the effective prevention and treatment of the disease.

This study investigated and analyzed transcriptomic data from nasal polyp tissues and normal mucosal tissues. The analysis indicated that the DEGs played pivotal roles in nucleocytoplasmic transport, regulation of DNA methylation, nucleus organization, and other processes. A recent meta-analysis incorporating 65 studies demonstrated that there was a notable distinction in overall methylation between the CRSwNP and control groups.²¹ Moreover, these DEGs were evidenced to have a part to play in various pathways, such as those associated with influenza A, dopaminergic synapse, and HIF-1 signaling pathway. Existing research indicates that CRSwNP is possibly induced by hypoxia, a factor that can contribute to hyperpermeability of airway epithelium and epithelial mesenchyme through upregulation of the HIF1 pathway.^22,23 The application of three machine learning methods to screen DEGs and subsequent comparison with the validation set GSE136825 revealed a high degree of consistency in gene expression patterns of BTBD10, CYFIP1, ERAP1, GIPC1, NYNRIN and PEX6. However, these genes have not yet been thoroughly studied in CRSwNP. Furthermore, an immuno-infiltration analysis was conducted on these genes. Notably, the results indicated reduced infiltration levels of type 2 T helper cells, activated CD4 T cells, memory B cells, mast cells, and natural killer T cells in the CRSwNP group, and increased infiltration levels of Immature dendritic cells and effector memory CD8 T cells. This finding is consistent with previous studies on NP.^9,24–27 Finally, following the successful establishment of a mouse mode, RT-qPCR validation was conducted on mouse nasal mucosal samples, which revealed that a higher degree of mRNA expression was noted for BTBD10, GIPC1, NYNRIN and PEX6 in the CRSwNP group than those in the control group. It is postulated that these genes may serve as potential new diagnostic or therapeutic targets.

The PEX6 gene is involved in the synthesis and maintenance of peroxisomes (peroxisome). Pex1 and Pex6 are ATPases related to a variety of cellular operations (AAA-ATPases) and are indispensable for the biogenesis and maintenance of peroxisomes. Being the exclusive Pex proteins capable of utilizing energy, Pex1 and Pex6 initiate the uptake of matrix enzymes into peroxisomes and proactively inhibit the deterioration of peroxisomes. Peroxisomes play a significant role in intracellular lipid metabolism, oxidative reactions, and lipid protonation processes.^28–31 In humans, mutations that result in the dysfunction of PEX1, PEX6, and PEX26 are the most widespread triggers for infrequent genetic conditions (known as peroxisome biogenesis disorders (PBD)).^31,32 Previous studies have indicated that aberrant regulation of PEX6 is a key contributor to diseases such as Zellweger syndrome,³³ Refsum disease,³⁴ and Usher syndrome.³⁵ In conclusion, abnormal expression of the PEX6 gene may result in dysfunctional peroxisomes, which may subsequently affect intracellular metabolic balance. The implications of such dysfunction are particularly relevant in the context of CRSwNP. The altered intracellular environment, characterized by a deranged metabolic state, can create conditions that are conducive to inflammation and tissue damage, which potentially stimulate the development and progression of CRSwNP. For instance, the accumulation of reactive oxygen species (ROS) due to impaired peroxisomal antioxidant defenses can lead to oxidative stress, a known contributor to the pathogenesis of nasal polyps.^36,37 Moreover, the dysregulation of lipid metabolism that might occur due to peroxisomal dysfunction could contribute to the hyperplasia of the nasal mucosa,^37,38 which might lead to a chronic inflammatory state that perpetuates the symptoms of CRSwNP. Therefore, the abnormal expression of the PEX6 gene and the subsequent impact on peroxisomal function could represent a significant genetic factor in the complex etiology of CRSwNP.

GIPC1, a protein situated in the cytoplasm, is found at the periphery of the membrane, and acts as an adaptor protein connecting receptor interactions with intracellular signaling pathways (including cell cycle regulation). It plays a role in the transduction of various signaling pathways such as the PDGFR/PI3K/AKT signaling pathway and is essential in the context of integrin re-uptake during cell migration, angiogenesis, and cytokinesis.^39–41 Abnormal expression of GIPC1 may disrupt the equilibrium of intracellular signaling pathways, accordingly influencing physiological mechanisms including cell differentiation, multiplication, and apoptotic cell death. Upregulation of GIPC1 in breast cancer, ovarian cancer, gastric cancer, and pancreatic cancer facilitates tumor proliferation and invasion, and inhibition of GIPC1 may holds potential for therapeutic intervention of human cancers.^39,41,42 The integration and regulation within the same tumor environment may be dependent on the intensity and time frame of interaction with GIPC.⁴³ Additionally, some research suggests that dysfunction of the GIPC1 protein may influence the pathogenesis of oculopharyngodistal myopathy.⁴⁴ In CRSwNP, GIPC1 may contribute to the disease process by participating in the imbalance of intracellular signaling pathways and affecting core physiological processes, encompassing cell differentiation, proliferation, and apoptosis. The dysregulation of these processes due to GIPC1’s involvement could lead to the inflammation, tissue hyperplasia, and structural changes observed in CRSwNP, thereby contributing to the chronicity and severity of the condition. Understanding the specific mechanisms by which GIPC1 influences these pathways could offer new insights into the pathobiology of CRSwNP and potentially lead to novel therapeutic strategies for its management.

BTBD10, designated as glucose metabolism-related protein 1 (GMRP1) as well, represents a fresh entrant of the KCTD family. BTBD10 harbors a BTB/POZ domain and has been identified to activate Akt by reducing dephosphorylation mediated by protein phosphatase 2A. It is instrumental in the processes of cell division and proliferation. The protein synthesized by this gene may participate in the division of the cell nucleus and the maintenance of chromosome stability, which is crucial for the growth and development of tissues and organs.^45–47 Earlier investigations have demonstrated that the overexpression of BTBD10 in human glioblastoma cells exerts antiproliferative effects and promotes apoptosis.⁴⁸ Additionally, BTBD10 overexpression has been revealed to facilitate the proliferation of pancreatic β cells.⁴⁹ Moreover, elevated BTBD10 expression in hepatocellular carcinoma tissues is associated with a worse prognosis.⁵⁰ In CRSwNP, the dysregulation of BTBD10, a protein pivotal to various cellular processes, has been observed to have significant implications on the behavior of epithelial cells. This abnormal expression of BTBD10 does not merely perturb the cells’ routine functions; it can fundamentally alter their ability to proliferate and to engage in the repair mechanisms necessary for maintaining tissue integrity. The consequences of such alterations are profound, as they lead to an exacerbated disease process in CRSwNP. With epithelial cells less capable of regenerating and repairing themselves, the cycle of inflammation and tissue damage is perpetuated, resulting in a chronic state that is epithelial barrier deficiency, which not only contributes to the persistence and severity of the polyps but also impacts the overall pathophysiology of the disease,^20,51 leading to a more complex clinical picture that challenges therapeutic intervention strategies.

NYNRIN (known for its inclusion of an NYN domain and retroviral integrase) is a common RNA-binding protein (RBP). Previous studies have demonstrated high expression of NYNRIN in pediatric acute myeloid leukemia with a high-risk profile and hereditary nephroblastoma. However, the function of NYNRIN remains largely unknown. Although the NYN domain is believed to serve a function in RNA processing events, and NYNRIN is associated with microRNA-mRNA regulation,^52,53 the functional mechanisms of this gene remain to be explored. Further investigation is warranted to ascertain the association between the NYNRIN gene and CRSwNP, with the objective of elucidating its role in the disease process.

It is well-established that hypoxia and immune mechanisms may be two key factors in the progression of CRSwNP, and they are inevitably intertwined. Currently, it is currently acknowledged that the involvement of hypoxia in the disease process of CRSwNP has been established, as oxygen content in inflammatory sinus cavities of CRSwNP has been found to be significantly reduced.⁵⁴ A recent study demonstrates that under hypoxic conditions, on the one hand, hypoxia regulates the multiplication and differentiation of macrophages, eosinophils, basophils, as well as mast cells in the sinus mucosa, thereby influencing the inflammatory status of CRSwNP via the adjustment of T cells and B cells.^54–57 On the other hand, the stable expression of hypoxia-inducible factors-1 (HIFs-1) α and HIFs-2α engages in the immune response and inflammatory pathways of CRSwNP. This affects the differentiation of nasal epithelial cells, aids in fibroblast proliferation, and triggers epithelial-mesenchymal transition (EMT) and tissue remodeling, which in turn mediates the formation of NP.^54,58–60 Furthermore, our ssGSEA results also indicate that type 2 T helper cells, activated CD4 T cells, memory B cells, and immature dendritic cells bear a close association with the pathogenesis of CRSwNP. However, the specific mechanisms require further experimental verification in future studies. Concurrently, the findings of our functional enrichment analysis indicate that DEGs are primarily enriched in the pathways related to the immune system and hypoxia. Consequently, therapy strategies that address the immune and hypoxic pathways are of importance and promise.

In this study, WGCNA analysis and three machine learning algorithms (LASSO, SVM-RFE, and RF) were employed to identify hub genes associated with CRSwNP, thereby reducing potential bias to the greatest extent. To validate the key genes, animal modeling experiments were conducted, which identified four hub genes as potential new targets for the diagnosis and treatment of CRSwNP.

It must be pointed out, however, that this study is not without its constraints. The sample size of the GSE136825 and GSE107624 validation cohorts was relatively limited, and thus, the results of the analyses should be validated in larger cohorts to determine the reproducibility of the study findings. While several hub genes that may be established as fresh targets for CRSwNP diagnosis and treatment were assessed using bioinformatics analysis, larger-scale prospective studies are required to validate our conclusions. Furthermore, more laboratory and biological studies on mechanisms both in test tubes and in living organisms should be conducted to expand upon and validate the results obtained by experiments.

In future research, we anticipate more methods to identify appropriate solutions to address these issues. We need to conduct functional studies to confirm the roles of these genes, further utilize animal models to investigate the underlying mechanisms, and carry out clinical research to explore the roles of the identified genes as therapeutic targets or diagnostic markers.

Conclusion

In conclusion, this study implemented three machine learning algorithms in combination with WGCNA analysis and animal model validation to identify BTBD10, GIPC1, NYNRIN, and PEX6 as potential hub genes that may be used as promising targets for the diagnosis and treatment of CRSwNP. These genes may be involved in the occurrence and development of CRSwNP through hypoxia and immune mechanisms. In the future field of CRSwNP, we must explore the functions of these genes through studies, use animal models to understand their mechanisms, and large-scale perform clinical trials to assess their potential as treatments or diagnostics.