The Gene Expression Omnibus (GEO) database is an extensive public resource for gene expression data across various species (21). We obtained the transcriptome data from the GEO database, which can be accessed at https://www.ncbi.nlm.nih.gov/geo/. For our study, the datasets must satisfy the following criteria: 1) The AR group should comprise at least 20 participants; 2) The research must focus on Homo sapiens; 3) The raw or processed data must be publicly available; 4) The research must involve blood samples from both individuals with AR and healthy individuals. The datasets used in this study consisted of the following: GSE75011, which included 15 healthy control (HC) samples and 25 samples from individuals with AR, and GSE50223, which included 21 HC samples and 21 AR samples. GSE75011 was chosen as the test set, while GSE50223 was selected as the validation set. The data obtained from these GEO datasets underwent preprocessing and normalization using the “affy” R package.

Using the limma package, we performed a differential expression analysis on the GSE75011 dataset, comparing the HC and AR groups. The DEGs were determined based on the criteria of a p-value < 0.05 and a |log-fold change (FC)| > 0.5.

To determine the most significant pathways between the HC and AR groups, we utilized the ClusterProfiler package to perform Gene Set Enrichment Analysis (GSEA). We obtained the “h.all.v7.4.symbols.gmt” subset from the Molecular Signatures Database to assess pertinent pathways. Statistical significance was defined as a p value below 0.05. Subsequently, the gene set variation (GSVA) method was employed to assess the inflammatory score. To accomplish this, we employed the GSVA package to calculate the inflammatory response score for each sample within the inflammatory response gene set. The resulting inflammatory response score was then visually represented using a boxplot. We obtained a total of 200 inflammatory response-related genes (IRGGs) from the MSigDB (inflammatory response gene set). To obtain the gene expression profile of DIRGGs, we intersected these IRGGs with the DEGs identified from the GSE75011 dataset. By utilizing the “pheatmap” R package, we generated a heat map that illustrates the expression patterns of DIRGGs. We used the ClusterProfiler package to perform functional enrichment analyses on DIRRGs. Our analysis involved the use of the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO). GO function annotation includes three categories: biological process (BP), cell component (CC), and molecular function (MF). The statistical significance was determined by enrichment, with a significance level set at p < 0.05. Ultimately, the results were visually represented using the ggplot2 package.

The STRING database (https://string-db.org/) is a valuable resource that allows users to investigate protein interactions (22). We rely on the STRING database to construct PPI networks, with a confidence score of 0.4 serving as the threshold for determining the reliability of these interactions. CytoHubba offers 11 different topological analysis techniques, encompassing Maximal Clique Centrality Degree (MCC), Density of Maximum Neighborhood Component, Maximum Neighborhood Component, Edge Percolated Component, and six centralities (Stress, Betweenness, Radiality, Closeness, EcCentricity, and Bottleneck). Out of these eleven methods, the recently introduced MCC method demonstrates superior precision in predicting essential proteins from the PPI network, setting it apart from the rest (23). We utilized the cytoHubba plugin within the Cytoscape software to identify hub genes using the MCC method. The top 10 genes with the highest MCC scores were chosen as hub genes, and these selected genes were further investigated in subsequent analyses (24).

To identify potential candidate genes for diagnosing AR, two machine learning algorithms were used. The randomForest package in R was used to construct a model and calculate the average error rate for all DIRRGs. Another random forest (RF) model was created to determine importance values by reducing accuracy. Genes with an importance value > 1 were identified as hub genes for future model development (25). The top 8 genes were chosen as promising candidates. The glmnet package was used for LASSO analysis, with penalty parameters for 10-fold cross-validation. This approach surpasses traditional regression analysis for high-dimensional data assessment (26). Subsequently, we obtained diagnostic genes by identifying the genes based on their scores in MCC, RF, and LASSO, and then intersecting the resulting gene lists.

The construction of a nomogram is extremely advantageous for the diagnosis of clinical AR. To create the nomogram, we employed the rms R package and considered the candidate genes (27). Each candidate gene was given a score known as “points”, while the “Total Points” represents the total score of all the genes. The diagnostic efficacy of both the candidate genes in AR diagnosis was evaluated by establishing a ROC curve. Additionally, the model’s accuracy was assessed using decision curve analysis (DCA) and calibration curves.

To conduct single gene GSEA analysis, we utilized the GSEA software to classify the samples into high and low expression groups, using the median values of diagnostic gene expression levels as the criteria. In order to explore the underlying molecular mechanisms associated with gene phenotypes, we obtained the “c2.cp.kegg.v7.4.symbols.gmt” subset from the Molecular Signatures Database. The minimum gene set was set to 5, while the maximum gene set was set to 5000. Additionally, one thousand resamplings were performed. A p value lower than 0.05 was considered to be statistically significant.

The ssGSEA algorithm utilized immune gene sets that incorporated genes associated with various immune cell checkpoints, pathways, functions, and types. To comprehensively evaluate the immunological attributes of each sample, we applied the ssGSEA algorithm using R packages such as limma, GSEABase and GSVA (28). Furthermore, the ggplot2 package was utilized to investigate the association between the expression of diagnostic genes and the infiltration of immune cells. The findings were visually presented in a lollipop chart.

The ConsensusClusterPlu package was employed in the study to conduct unsupervised cluster analysis using the expression profile of 22 DIRGGs. This facilitated the categorization of AR patients into distinct subgroups and the identification of the most suitable number of clusters.

The Connectivity Map (CMap) is a database of expression profiles that harnesses cellular responses to disturbances in order to uncover potential functional connections between treatments, genes, and diseases (29). In this study, we performed CMap analysis to anticipate small-molecule compounds that target subtypes associated with inflammatory response. We extracted drug signatures from the Connectivity Map database (CMap, https://clue.io/). Our input data consisted of the top 100 up-regulated and 100 down-regulated genes in the two subtypes. The final analysis result is assigned a score ranging from -100 to 100. In the case of small-molecule drugs, a negative score suggests the ability to reverse gene expression, which signifies their potential therapeutic significance. We chose the top five small-molecule compounds based on their lowest CMap score.

The binding affinity between signature gene proteins and small-molecule drugs was investigated by molecular docking technique (30). The protein structure of signature genes was obtained from the PDB database, while the molecular structures of small-molecule drugs were acquired from PubChem. Ligands and receptors were then processed using AutoDockTools, primarily through the addition of hydrogen atoms and identification of active pockets. The binding modes between the candidate protein and small-molecule drugs were subsequently analyzed using AutoDock Vina software.

A total of 16 individuals participated in the study, with an equal distribution of 8 HC and 8 patients diagnosed with AR. The blood samples were collected from the Affiliated Huai’an Hospital of Xuzhou Medical University. Prior to the collection, all participants provided written informed consent, and the study protocol received approval from the Ethics Committee of the Affiliated Huai’an Hospital of Xuzhou Medical University (ethical approval number: 2022029). We collected peripheral blood in tubes containing citrate, which were then stored at 4°C until ready for use.

RNA was extracted from blood samples using the TRIzol reagent (Invitrogen, CA, USA) following the manufacturer’s instructions. Then, the RNA underwent spectrophotometric quantification at a wavelength of 260 nm, and its purity was assessed by determining the absorbance ratio at 260/280 nm. The absorbance ratios (260/280 nm) of the RNA samples are fall within the range of 1.8 to 2.1. 2 μg of RNA was then subjected to reverse transcription using the First Strand cDNA Synthesis Kit (Invitrogen, CA, USA). qRT-PCR analysis was performed using the SYBR qPCR Master Mix (Sigma, MO, USA). The parameter settings during amplification are as follows: initial denaturation at 95°C for 10 minutes, followed by 40 cycles consisting of denaturation at 95°C for 10 seconds, annealing at 60°C for 30 seconds, extension at 72°C for 1 second, and final cooling at 40°C for 30 seconds. The expression levels of the diagnostic signature genes were measured using the Roche LC480 Real-Time PCR System (Roche). The internal control for mRNA was β-actin. The primers were presented in Supplementary Table S1.

The flowchart illustrating the analysis process for this study is presented in Figure 1. A comparative analysis was conducted between samples from HC and AR, resulting in the identification of 1323 DEGs. Among these DEGs, 307 were upregulated, whereas 1016 were found to be downregulated (Figure 2A). The heat map displayed the level of DEGs in the top 50 between the two groups. The majority of these genes exhibit a downward trend (Figure 2B). The GSEA analysis revealed a significant enrichment of the inflammatory response gene set in the HC group (Figure 2C). Moreover, the GSVA algorithm was utilized to calculate the inflammatory response score, highlighting a significantly lower score in the AR group compared to the HC group (Figure 2D). Based on our findings, it becomes evident that the development of AR is closely associated with the inflammatory response. Therefore, our study aims to delve deeper into the role of IRGGs in AR, shedding light on their significance in this process.

Enrichment analysis was performed on the overlapping genes between DEGs and IRRGs. In this study, we identified a total of 22 DIRRGs, as shown in Figure 3A. Our findings demonstrated that the majority of these DIRRGs, such as PLAUR, NFKBIA, RGS16, GPR132, PTGIR, HIF1A, CXCR6, CDKN1A, LIF, CD55, PDE4B, GNA15, PTGER4, NAMPT, NLRP3, CCRL2, ICOSLG, and CCL20, exhibited a significant decrease in samples from AR patients compared to those from HC, as illustrated in Supplementary Figure S1. According to the data presented in Figures 3B, C, the results of the GO analysis revealed that these genes were significantly associated with BP, including response to lipopolysaccharide, cellular response to chemokine, and chemokine-mediated signaling pathway. Additionally, in terms of CC ontology, these genes were found to be involved in the inflammasome complex. Moreover, the MF analysis identified C-C chemokine binding, C-C chemokine receptor activity, and icosanoid receptor activity as significant terms for these genes. Furthermore, the KEGG analysis demonstrated that these genes were associated with pathways such as NOD-like receptor signaling pathway, chemokine signaling pathway, TNF signaling pathway, etc. We have also identified key players involved in the development of chronic myeloid leukaemia - CDKN1A, MYC and NFKBIA. These genes also intersect with key signaling pathways that control inflammatory responses and cell communication. In particular, the genes MEFV, NAMPT, NFKBIA and NLRP3 were involved in the NOD-like receptor signaling pathway, which orchestrates the body’s innate immunity. Meanwhile, CCL20, LIF and NFKBIA emerged as components of the TNF signalling pathway, a key mediator of inflammation and apoptosis. Similarly, the JAK-STAT pathway, a conduit for a variety of cytokines and growth factors, was influenced by the activities of CDKN1A, LIF and MYC. Finally, the chemokine signaling pathway, a key player in directing immune cell traffic, involves CCL20, CXCR6 and NFKBIA. Each of these pathways, shown in Figure 3C, reveals a complex tapestry of genetic interactions that are critical to our understanding of the mechanisms of AR.

A PPI network was established, illustrating the interaction between 19 different genes (Figure 4A). The ranking of these genes based on their MCC score can be observed in Figure 4B. To identify signature genes in patients with AR, two machine algorithms were employed. The random forest analysis successfully identified 9 signature genes with a relative importance exceeding 1, as demonstrated in Figure 4C. Additionally, the LASSO analysis specifically selected 8 signature genes, as illustrated in Figure 4D. Figure 4E displayed a Venn diagram depicting the overlap between the top 10 genes identified using MCC, the 8 significant genes identified using RF, and the 8 potential candidate genes identified using LASSO. For analysis and validation purposes, four signature genes (NFKBIA, HIF1A, MYC, CCRL2) were identified from this intersection.

In the AR group, the levels of NFKBIA, HIF1A, and CCRL2 expression were lower compared to the HC group. Conversely, the mRNA expression of MYC was higher in the AR group than in the HC group (Figure 5A). In order to improve the accuracy of predicting the development of AR patients, a nomogram has been created which incorporates the analysis of four specific genes (Figure 5B). By analyzing the receiver operating characteristic (ROC) curve, it was observed that the model performed exceptionally well, with a significant area under the curve (AUC) value of 0.931 (Figure 5C). The results from the calibration curve (Figure 5D) further validated the remarkable precision of the nomogram model in predicting outcomes for patients with AR. Moreover, the decision curve analysis demonstrated the potential advantages of employing the nomogram model in patients with AR, as shown in Figure 5E. Moreover, the model’s reliability was assessed by using an external validation dataset (GSE50223) (Supplementary Figure S2). The results demonstrated that the nomogram achieved a remarkable diagnostic accuracy for AR, highlighting its efficacy.

To further validate the diagnostic value of the four markers in diagnosing AR, we gathered clinical blood samples. With these samples, we aimed to confirm the effectiveness of these markers in accurately identifying AR cases. In Figures 6A–D, the AR group had lower expression of NFKBIA, HIF1A, and CCRL2 compared to HC, while the AR group had higher expression of MYC (p < 0.05 or p < 0.01 or p < 0.001). Based on our research, we discovered that the above-mentioned genes associated with the inflammatory response exhibit a strong diagnostic capability, making them potential biomarkers for diagnosing AR.

We conducted a detailed examination of potential signaling pathways linked to signature genes using GSEA. The results obtained from single gene GSEA indicated that the NFKBIA high-expressed phenotype was enriched with MAPK signaling pathway, NOD like receptor signaling pathway, B cell receptor signaling pathway, chemokine signaling pathway, endocytosis, and apoptosis (Figure 7A); the HIF1A high-expressed phenotype was enriched with T cell receptor signaling pathway, B receptor signaling pathway, chemokine signaling pathway, chronic myeloid leukemia, TGF beta signaling pathway (Figure 7B); the MYC low-expressed phenotype was enriched with chemokine signaling pathway, toll like receptor signaling pathway, and MAPK signaling pathway (Figure 7C); the CCRL2 low-expressed phenotype was enriched with alpha linolenic acid metabolism, primary bile acid biosynthesis, linoleic acid metabolism, and sulfur metabolism (Figure 7D). These findings provided additional evidence that the NFKBIA, HIF1A, and MYC genes are all intricately connected to the immune system and the body’s response to inflammation.

In addition, we employed the ssGSEA algorithm to quantify the enrichment scores of immune cells in the blood of AR. The results showed a notable disparity in immune cell infiltration between the HC and AR groups, with the HC group exhibiting significantly higher levels of Treg and TFH cells (Figure 8A). By analyzing the correlation between immune cells, we have identified certain genes that may be involved in the progression of AR. These genes exert their influence by regulating immune cells like CD8 T cells, mast cells, neutrophils, B cells, and THF (Figures 8B-E). In conclusion, our study provided valuable insights into the underlying mechanisms of AR and the role of immune cell infiltration in its progression. The identification of signature genes associated with specific immune cells not only enhances our understanding of the AR but also opens up new avenues for targeted therapeutic interventions. Further research is needed to validate these findings and explore the potential clinical applications in the management of AR.

We also classified AR patients into different inflammatory phenotypes using unsupervised cluster analysis. This allows us to better understand the variation and characteristics within the patient population. Based on the expression profiling of 22 DIRGGs, a consensus clustering approach was used to cluster AR patients. The results revealed that the optimal number of subtypes was determined to be 2 (Figure 9A). Notably, significant heterogeneity in gene expression was observed between these two subtypes. A comparative analysis of samples from subgroups C1 and C2 found a total of 1463 differentially expressed genes. Among these, 463 genes were observed to be upregulated, while 1000 genes were observed to be downregulated (Figures 9B, C). Furthermore, the results of KEGG and GO analysis indicated that these differentially expressed genes were predominantly enriched in pathways related to immune and inflammatory processes. These pathways encompassed autophagy, infection by human T cell leukemia virus 1, chronic myeloid leukemia, proteolysis, cytokine production involved in immune response, myeloid cell development, T helper 2 cell differentiation, and other related pathways (Figures 9D, E).

We also performed GSEA analysis on the entire set of genes linked to two distinct inflammatory patterns. Our results demonstrated notable variations in certain pathways, particularly in the TGF beta signaling, TNFA signaling via NFKB, P53 pathway, etc (Figure 10A). In addition, GSVA results revealed that growth hormone receptor signaling via JAK-STAT, negative regulation of autophagy, acute myeloid leukemia, cell growth, TGF beta signaling, positive regulation of inflammatory response to antigenic stimulus and P53 pathway were activated in the C2 subgroup, which were consistent with the results of KEGG and GSEA (Figures 10B, C). Collectively, our findings from both GSEA and GSVA analysis, along with the functional enrichment analysis, provide comprehensive evidence supporting the activation of specific pathways in the C2 subgroup. These findings shed light on the molecular mechanisms underlying the inflammatory patterns observed in this subgroup and may pave the way for the development of targeted therapeutic approaches for AR associated with these activated pathways.

Through the utilization of CMap analysis, we aim to uncover potential therapeutic drugs that target AR inflammation-related subtypes. A small molecule compound with a lower CMap score is more likely to possess the capability to treat the AR. As shown in Supplementary Table S2, Levetiracetam, oxybutynin, metyrapone, idebenone, and LM-1685 emerged as the five small-molecule drugs with the most favorable CMap score for the C1 subtype (Figure 11A). Conversely, triptolide, daunorubicin, dactinomycin, tipifarnib, and atorvastatin ranked as the top five small-molecule drugs with the lowest CMap score for the C2 subtype (Figure 12A).

For the molecular docking, we selected five drug candidates and identified the signature genes as NFKBIA, HIF1A, MYC, and CCRL2. The molecular docking results were presented in Figures 11B, 12B. Among these small molecules, LM-1685 and dactinomycin exhibited the strongest binding affinity to the target proteins. Subsequently, we utilized PyMOL software to visualize the molecular docking results, as shown in Figures 11C–F, 12C–F.