The flowchart depicted in Figure 1 outlines the systematic approach undertaken in this study. All datasets employed for this research are accessible to the public, ensuring transparency and reproducibility. The transcriptome data were sourced from the publicly accessible Gene Expression Omnibus (GEO) database at https://www.ncbi.nlm.nih.gov/geo/. High-throughput sequencing data from GSE176205 and GSE167199 were corrected for batch effects using the “sva” package(version 3.42.0) in R (23). The sequencing platform for GSE176205 was GPL20301 Illumina HiSeq 4000 (Homo sapiens), which included 6 IDD and 3 standard control NP samples. The sequencing platform for GSE167199 was GPL24676 Illumina NovaSeq 6000 (Homo sapiens), including 3 IDD and 3 standard control NP samples. GSE167931 was utilized as an additional dataset for validation purposes, sequenced on the GPL20795 HiSeq X Ten (Homo sapiens) platform, and included 5 IDD and 4 standard control NP samples.

A single-cell sequencing dataset was retrieved from the GEO database, which contained 16 samples (8 IDD and 8 standard control NP samples). The “Seurat” R package(version 4.2.0) was employed in the R to import the raw data within GSE153066 from the database (24). Filtering out low-quality cells and genes involved applying stringent criteria to ensure data quality: 1) Cells that failed to express at least 200 genes were excluded. 2) Cells with a mitochondrial gene expression percentage below 25% were also retained. 3) Cells with fewer than 4,000 expressed genes were retained. 4) Cells with a UMI reading of fewer than 10,000 were retained. The data were normalized using the “Seurat” R package. Following normalization, balancing average expression and dispersion was employed to pinpoint highly variable genes. Principal component analysis (PCA) yielded significant PCs, which served as input for graph-based clustering. A harmony approach was used to eliminate batch effects across samples. For clustering, the “FindClusters” function, leveraging the optimized Shared Nearest Neighbor (SNN) modular clustering algorithm, was applied to identify 16 distinct clusters with a resolution of 0.5. The “RunUMAP” function was used for Uniform Manifold Approximation and Projection (UMAP). UMAP-1 and UMAP-2 visualized cell aggregation, subsequently enabling the identification of cell clusters based on their type-specific biomarkers. The “RunTSNE” function was then used for t-distributed stochastic neighbour embedding(t-SNE). t-SNE-1 and t-SNE-2 visualized T cell aggregation.

ssGSEA represents an extension of Gene Set Enrichment Analysis(GSEA) and it was employed to calculate enrichment scores for each cell type in conventional transcriptome data based on a list of marker genes (25). Based on a list of marker genes for each cell type, we used “GSVA” R package (version 1.42.0) to perform ssGSEA analysis to calculate the enrichment fraction of each cell type in transcriptome data.

Support Vector Machine-Recursive Feature Elimination (SVM-RFE), a machine learning approach, trains subsets of features across categories to refine the set and pinpoint the most predictive genes. Utilizing the “glmnet” package in R (version 4.1.4), LASSO regression was implemented to compute and select linear models, effectively retaining only informative variables. Random forest analysis was conducted using the “RandomForest” function, where the minimum error was chosen as the mtry node value, and the stabilized image value was selected as the ntree. The results of SVM-RFE, LASSO regression, and Random forest were combined to select the key immune cells in this study.

Immune Response Enrichment Analysis (IREA) is a method Ang Cui et al. proposed to infer cytokine activity and immune cell polarisation during immunization (26). In a cell type-centric view, IREA uncovers over 66 distinct cytokine-mediated polarisation states within immune cell types, encompassing novel and previously undescribed states. IREA was used to calculate the cytokines that primarily generate a response based on the differentially expressed genes (DEGs) of vital immune cells and the polarisation states of critical immune cells in IDD.

Using the “CellChat” (version 1.1.3) R package, CellChat objects were constructed (27). Cell-cell communication was analyzed using the “CellChatDB.human” ligand-receptor interaction database as the reference. The “mergeCellChat” function was employed to combine the CellChat objects of each group, enabling a comparison of the interaction number and the interaction strength. The “netVisual_diffInteraction” function was used to visualize differences in the number or strength of interactions between groups and different cell types. Finally, the “netVisual_aggregate” function was used for visualization.

The DEGs were identified utilizing the “limma” R package (version 3.50.0) (28). The screening criteria were |log2Fold Change|>1 and adjusted p-value <0.05. Heatmaps were created employing the “pheatmap” R package(version 1.0.12) and clustered using euclidean distance and hierarchical clustering methods. Further, we explored the differences between the two groups in single-cell data. The “FindAllMarkers” function in the “Seurat” R package was employed to identify DEGs of vital immune cells between two groups.

A PPI network was constructed by utilizing the Search Tool for the Retrieval of Interacting Genes (STRING) online resource (29). Sub-networks were constructed for proteins with interaction scores greater than 700, from which proteins with direct interactions with intersecting DEGs were extracted and analyzed. Afterward, PPI network analysis was executed using the Cytohubba plug-in (30). Ten algorithms (Betweenness, BottleNeck, Closeness, Degree, EcCentricity, EPC, MCC, MNC, Radiality and Stress) were used to screen the top 100 proteins in terms of importance. The intersecting proteins from these algorithms and the genes identified in the intersecting DEGs were designated as key genes.

GO enrichment analysis includes three domains: biological process (BP), molecular function (MF), and cellular component (CC) analysis (31). The KEGG is a bioinformatics resource for identifying enriched and significantly altered metabolic pathways within gene lists (32). The genes of interest in IDD were analyzed for GO and KEGG enrichment using the R package “clusterProfiler” (version 4.2.2) (33).

The GeneMANIA website (http://genemania.org) forecasts associations between functionally analogous genes and pivotal genes (34). Gene interaction networks for key genes were constructed using the GeneMANIA website.

To explore the disparities in biological functionalities between the two groups, the gene set “h.all.v2023.1.Hs.symbols.gmt” sourced from the MSigDB was employed as the reference gene set. GSVA was performed using the R package “GSVA” (version 1.42.0).

Since the mechanism of action of competing endogenous RNA (ceRNA) in IDD remains unclear, we used miRTarBase (https://mirtarbase.cuhk.edu.cn/~miRTarBase/miRTarBase_2022/php/index.php), starbase3.0 (https://rnasysu.com/encori) and miRDB databases (https://mirdb.org/index.html) to back-predict microRNAs of key genes and to predict lncRNAs of common microRNAs of key genes, constructing a ceRNA network (35, 36).

The open-source platform (https://starbase.sysu.edu.cn/tutorialAPI.php#RBPTarget) was used to investigate mRNA and RNA-binding protein (RBP) expression associations.

The mRNA-TF interaction relationships between key genes were constructed using the TRRUST(Transcriptional Regulatory Relationships Unraveled by Sentence-based Text Minmining) database. Networks were constructed using Cytoscape(version 3.9.1).

We acquired human NPCs from Procell Life Science & Technology Co., Ltd (CP-H097, Wuhan, China). We cultured in F12/DMEM medium (Gibco, USA) supplemented with 1% Penicillin-Streptomycin Solution (NCM Biotech, China) and 10% fetal bovine serum (FBS; Gibco, USA). As described in previous literature, NPCs were induced using 100 µM TBHP (Sigma-Aldrich, St. Louis, MO, USA) for 4 hours to construct the IDD model in vitro (37). Cells were transfected with small interfering RNA (si-RNA) to knock down IRF7 (ENSG00000185507, Gene ID: 3665) using the riboFECT CP Transfection Kit (RiboBio, Cat. No C10511-05, China).

Total RNA was extracted from the cells using TRIzol (Thermo Fisher, USA) and chloroform (Thermo Fisher, USA). The purity and concentration of RNA samples were quantitatively assessed utilizing a Nanodrop One spectrophotometer (Thermo Fisher, USA). cDNA synthesis from the extracted RNA was performed using the PrimeScript RT reagent. Reverse transcription reagents were purchased from TaKaRa, and an RT-qPCR assay was conducted using an ABI 7500 Real-Time PCR machine (Thermo Fisher, USA). Supplementary Table 1 provides the primer sequences employed in the study.

The Ethics Committee of Nanchang University approved all experimental procedures involving animals, and all procedures followed the ARRIVE guidelines. Male Sprague-Dawley (SD) rats, aged 12 weeks, were sourced from SpePharm Biotechnology Ltd. (Beijing, China). A rat tail puncture model was established according to the protocol described in the literature (38). Rats were administered 4% sodium pentobarbital (1mL/1kg) to induce anaesthesia. Following skin sterilization, a 21-gauge needle was employed to make a vertical puncture into the Co4/5 disc of the tail, reaching a depth of approximately 5 millimeters. Upon completing the puncture, the needle was rotated in a full 360-degree motion and maintained in position for one minute. On postoperative days 1, 7, and 14, si-RNA and control agent si-NC 5 nmol (10 µL) were delivered directly into the disc along the original puncture channel. The health status of the rats was monitored daily after surgery. The flowchart of the animal experiment was performed by Figdraw (www.figdraw.com).

X-rays and MRIs of the caudal intervertebral discs of rats were performed at 4 weeks postoperatively. The disc height index (DHI) score was utilized to quantify the structural integrity of the discs. This score was derived from the ratio of the intervertebral space height to the height of the neighbouring vertebral bodies observed in X-ray images. Additionally, the severity of IDD was assessed using MRI images and analyzed with the Pfirrmann grading system as a standardized reference.

4 weeks post-surgery, the rats were euthanized, and their disc samples were fixed and immersed in 4% paraformaldehyde. Following the EDTA decalcification solution, the disc samples were dehydrated, embedded in paraffin, and sectioned serially. NP morphology was examined using hematoxylin-eosin (HE) and Safranin O-Fast Green (SO&FG) staining.

Sections were treated with a 3% hydrogen peroxide solution to block endogenous peroxidase activity, followed by 3% bovine serum albumin (BSA; Servicebio, Cat no: GC305010). Primary antibodies used were COL2A1(Proteintech,Cat No.28459-1-AP), MMP13(Proteintech,Cat No.18165-1-AP), NLRP3(Proteintech,Cat No. 19771-1-AP) and IL-1β(Proteintech,Cat No. 16806-1-AP). These were incubated overnight at a temperature of 4°C. HRP-coupled secondary antibodies (Servicebio, Cat no: GB23303) were then used for the corresponding primary antibodies. Finally, the samples were stained using diaminobenzidine and hematoxylin to label the nucleus.

The NPCs were fixed with 4% Paraformaldehyde Fix Solution (Biosharp, China, Cat no: BL539A) for 15 minutes. The cells were blocked with 3% BSA at room temperature for 30 minutes. Then, the cells were incubated overnight with antibodies against COL2A1, MMP13,NLRP3 and IL-1β. After washing the samples, the secondary antibody was added and incubated at room temperature for 50 minutes (Servicebio, Cat no:GB21303). The samples were rewashed, and DAPI staining solution was added for incubation at room temperature in the dark for 10 minutes.

Based on the top 20 up-regulated and 20 down-regulated genes with significant differences between the IDD group and the control group, small molecule compounds that might be effective in treating IDD were predicted using the Connectivity Map (CMAP) database. Eight compounds with the highest scores as potential drugs were selected, with toxic compounds excluded. We downloaded SDF format files for potential drugs from PubChem, PDB files of target proteins in the PDB database, and dehydrated and deliganded the proteins in PyMOL. Proteins were hydrogenated in AutoDockTools (version 1.5.7). Molecular docking was then performed in AutoDockTools. The strength of the binding energy serves as an indicator of the probability of the receptor and ligand binding together, with a decrease in binding energy corresponding to an increase in the affinity. The lower the binding energy, the more stable the conformation of the receptor became. The results of molecular docking were visualized in PyMOL.

Statistical analyses were performed utilizing R software (version 4.1.2). Spearman’s correlation test were employed to examine associations between two variables. The Wilcoxon rank-sum test was applied for comparisons between two groups, while the Kruskal-Wallis test was utilized for comparisons involving three or more groups. A statistical significance threshold was set at a p-value of less than 0.05.

After initial quality control and doublet removal of IDD and control samples, 65,165 cells were obtained from the single-cell dataset. All cells were aggregated into 16 clusters ( Figure 2A ). And annotated by cell-specific biomarkers, we finally found 9 cell types, namely: nucleus pulposus (NP, cluster 0,1,4,5,14), progenitor nucleus pulposus (Pro_NP, cluster 2, 11), Monocyte (cluster 3,6,15), Macrophage (cluster 7), Neutrophil (cluster 8), T cell (cluster 9), Erythrocyte (cluster 10), B cell (cluster 12), and Endothelial cell (cluster 13) ( Figures 2B, C ). The distribution and number of different cell types varied significantly between the two groups ( Figure 2D ). NPCs were most reduced considerably in the IDD group, while Pro_NP, Macrophage, Erythrocyte, B cell, and Endothelial cell were also reduced to varying degrees. Monocytes increased most significantly in the IDD group, and Neutrophil also increased in IDD. T cells were re-clustered and categorized into CD4⁺ T cells and CD8⁺ T cells based on their distinct expression patterns ( Figures 2E–H ).

The ssGSEA algorithm was used to calculate the enrichment scores of each cell type ( Supplementary Table 2 ). There was a notable decrease in the number of macrophages in the IDD group compared to the control group. In contrast, the IDD group’s monocytes count was significantly elevated compared to the control group. This is consistent with trends observed in the single-cell data ( Figure 3A ). The Random Forest, LASSO regression, and SVM-RFE were used to screen for critical immune cells. We selected the top 3 immune cells using the random forest algorithm based on Mean decrease accuracy (MDA ) and Mean decrease Gini (MDG), which were macrophages, monocytes, and neutrophils ( Figure 3B ). By LASSO regression analysis, we obtained the same three types of immune cells ( Figures 3C, D ). However, through the SVM-RFE method, we only screened for monocytes ( Figure 3E ). Finally, the immune cells detected by each method were intersected, concluding that the most critical immune cell in IDD was the monocytes ( Figure 3F ). The ROC curve showed that the monocytes enrichment score had good efficacy for distinguishing IDD samples from normal control samples (AUC=0.8333, Figure 3G ).

In order to clarify the changes of cytokines in monocytes during IDD and to explore the polarization state of monocytes after receiving different cytokines, we performed a comparative analysis of monocytes between different groups to explore the immune status of monocytes in IDD ( Figure 4A ; Supplementary Table 3 ). Subsequently, IREA was conducted based on significantly up-regulated differentially expressed genes in IDD. The cytokine enrichment plot showed that monocytes in IDD were predominantly enriched for cytokines such as IFN-a1 and IFN-b ( Figure 4B ; Supplementary Table 4 ). The cellular polarisation radar plot indicated that monocytes were primarily in the Mono-a polarisation state (type I interferon-induced polarisation), which may exacerbate the IDD’s inflammatory response ( Figure 4C ; Supplementary Table 5 ). When IDD occurs, IFN-a1 and IFN-b are significantly enriched in the intervertebral disc. And the monocytes are in a Mono-a polarized state, which may lead to the aggregation of inflammatory factors and matrix degradation in the intervertebral disc.

To investigate the cellular interaction network in IDD, the R package “Cellchat” was performed. Compared to the control group, the number and the strength of interactions between different cell types were elevated in the IDD group ( Figure 5A ). Besides, we showed the relationship of the number of interactions between different cell types in the IDD and the control group. Relatively little difference was found between the groups ( Figure 5B ). Then, we focused on monocytes and found that the biggest change occurred between monocytes and macrophages. And the communication strength between the two cells decreased greatly in IDD, indicating that the interaction between macrophages and monocytes is an important factor in the occurrence of IDD ( Figure 5C ).

Next, we visualized the relationship between immune cells and signaling pathways in the control and IDD group in the heatmap. It demonstrated that in the control group, ANNEXIN is the most important signaling pathway, but its communication signal decreases significantly after the occurrence of IDD ( Figure 6A ). To explore the reasons for this trend, violin plots were used to show the expression of the ANNEXIN pathway between the control and IDD groups. Compared with the control group, the ligand ANXA1 showed a decreasing trend in all cell types. And monocytes, with the highest expression of receptor FPR1, may be affected ( Figure 6B ). In addition, we also observed a significant decline in the CXCL signaling pathway in monocytes after the occurrence of IDD( Figure 6A ). Thus, the CXCL pathway-mediated communication between monocytes and other cells was mapped. The results showed that the CXCL pathway of monocytes mainly occurred in communication with macrophages and was significantly reduced in IDD ( Figure 6C ). Therefore, the reduction of CXCL signaling in monocytes and macrophages may exacerbate the progression of IDD.

2,500 DEGs were identified between IDD and control group, including 2,413 up-regulated genes and 87 down-regulated genes( Figures 7A, B ; Supplementary Table 6 ). By intersecting DEGs from the conventional transcriptome with DEGs from monocytes, 80 intersecting DEGs were identified as key genes ( Figure 7C ; Supplementary Table 7 ). GO enrichment analysis indicated that these genes play significant roles in BPs, such as regulation of the viral life cycle, regulation of viral processes, and viral genome replication. Additionally, CC, like tertiary granule lumen, focal adhesion, and cell-substrate junction, were enriched, as well as MF, such as the structural constituent of the ribosome ( Figure 7D ). KEGG analysis showed they were significantly enriched in Ribosome and Coronavirus disease - COVID-19 ( Figure 7E ; Supplementary Table 8 ). This suggests that ribosomes may be indirectly involved in the regulation of degenerative inflammatory responses by supporting the protein synthesis function of monocytes. The inflammatory environment of IDD activates monocytes, which synthesize and secrete a variety of inflammatory factors and degrading enzymes with the help of ribosomes, accelerating the destruction and degeneration of the intervertebral disc. Enrichment in viral biology suggests a role for the activation of inflammatory stress pathways in viral infection-mediated monocytes in IDD progression. Certain viruses, such as herpes simplex virus, can infect disc cells, causing direct cytopathy, inflammation, or apoptosis, which can accelerate disc degeneration (39). Coronavirus may also have a similar effect on disc degeneration.

A PPI network was created for the intersecting DEGs using the STRING online database. The top 100 proteins in terms of importance were screened using 10 algorithms, including Betweenness, BottleNeck, Closeness, Degree, EcCentricity, EPC, MCC, MNC, Radiality, and Stress. Six genes were identified through the intersection of these algorithms: POLR2A, RPLP0, JUN, CAT, IRF7, and RPS3 ( Figure 8A ). POLR2A, CAT, and IRF7 were included in the list of intersected DEGs, and RPLP0, JUN, and RPS3 were the genes that had a direct reciprocal relationship with the intersected DEGs. POLR2A, CAT, and IRF7 were selected as key genes for subsequent analysis. These essential genes showed up-regulation in the IDD group ( Figure 8B ). The ROC curves demonstrated POLR2A, CAT, and IRF7 all showed good efficacy for distinguishing IDD and standard samples ( Figures 8C–E ). Furthermore, the expression levels of these pivotal genes were confirmed in GSE167931. IRF7 were also significantly up-regulated. The expression of CAT in the IDD group is slightly upregulated, consistent with the trend in our analysis, which also shows potential as an IDD biomarker ( Figure 8F ). IRF7 showed promising efficacy for IDD and standard samples (AUC=0.9, Figure 8G ).

GSVA showed the differences between the IDD and control groups across 50 Hallmark signaling pathways ( Supplementary Table 9 ). In the IDD group, 3 pathways showed significant up-regulation and 15 pathways showed significant down-regulation ( Figure 9A ). IRF7 exhibits a significant negative correlation with HALLMARK_KRAS_SIGNALING_UP_PATHWAY, and down-regulation of HALLMARK_KRAS_SIGNALING_UP_PATHWAY has also been observed in IDD, suggesting that IRF7’s promotion of IDD may be mediated by inhibition of the KRAS pathway ( Figure 9B ). The activation of KRAS can upregulate the expression of the anti-apoptotic protein Bcl-2, thereby inhibiting apoptosis (40). In IDD, KRAS may be inhibited, leading to increased apoptosis of NPCs.

To delve deeper into the function of the key genes, a co-expression network was constructed utilizing the GeneMANIA database ( Figure 10A ; Supplementary Table 10 ). The correlation heatmap between hub genes is shown in Figure 10B , demonstrating that most key genes were closely correlated with each other. GO and KEGG enrichment analysis found that it is related to lipid metabolism and peroxide pathways. This suggests that three hub genes and closely interacting genes may be involved in IDD progression by mediating these processes ( Figures 10C, D ; Supplementary Table 11 ). Specifically, lipid metabolism disorders can lead to abnormal deposition of lipids in the intervertebral discs, triggering an inflammatory response (41). Lipid peroxidation accelerates disc cell damage and matrix degradation by increasing oxidative stress and promoting an inflammatory response (42). These processes interact with each other to form a vicious circle that ultimately leads to the progression of disc degeneration.

In order to further study the way of hub gene regulation, we constructed mRNA, TF, and ceRNA regulatory networks. An mRNA-miRNA-lncRNA interaction network was constructed using the key genes. Two miRNAs (miR-30a-5p and miR-30b-5p) were found to bind to CAT. A total of 23 lncRNAs were identified as target lncRNAs. The ceRNA network is shown in Figure 11A ( Supplementary Table 12 ). The corresponding mRNA/RBP pairs of hub mRNAs with pairwise information were searched and downloaded using the StarBase online database. Based on the relationships identified among the target genes from the online dataset, we constructed an intricate RBP-mRNA network ( Figure 11B ; Supplementary Table 13 ). Transcription factors (TFs) that bind to co-expressed genes were searched using the TRRUST database, identifying interaction relationship data for 2 target genes and 17 TFs ( Figure 11C ; Supplementary Table 14 ).

Using RT-qPCR analysis, it was noted that IRF7, POLR2A, and CAT exhibited significant upregulation in the TBHP-induced degeneration group, with IRF7 showing the most significant increase in expression ( Supplementary Figure 1 ). IDD can cause an imbalance between catabolism and anabolism in the ECM, leading to structural changes in the intervertebral disk (43). To investigate the potential of knocking down IRF7 to reverse IDD in vitro, we set up four groups of experiments: standard group, IDD group (TBHP), Si-NC group (TBHP + Si-NC), and Si-IRF7 group (TBHP + Si-IRF7). The efficiency of the IRF7 knockdown in NPCs was verified by RT-qPCR ( Supplementary Figure 1 ). It was found that knocking down IRF7 increased the expression of COL2A1 in the Si-IRF7 group while decreasing the expression of MMP13 compared to the Si-NC group ( Figures 12A, B, E ). It indicates that the knockdown of IRF7 can reverse the metabolic imbalance of the extracellular matrix (ECM). Besides,the knockdown of IRF7 decreased the expression of NLRP3 and IL-1β compared to the Si-NC group ( Figures 12C–E ). This indicates that knocking down IRF7 can alleviate the inflammation of IDD. In conclusion, these findings demonstrate that the knockdown of IRF7 can slow the progression of NPC degeneration.

Finally, we investigated the potential of knocking down IRF7 to reverse IDD in vivo. The flowchart of the animal experiment is shown in Figure 13A . Immunohistochemistry showed higher IRF7 expression in the Si-NC group than in the control group, and IRF7 was successfully knocked down in vivo ( Supplementary Figure 2 ). 4 weeks after the knockdown of IRF7, X-ray analysis showed an increased intervertebral space height in the Si-IRF7 group compared to the Si-NC group, with corresponding increases in DHI scores ( Figures 13B, C ). MRI results indicated significantly higher signal intensity in the Si-IRF7 group than in the Si-NC group, consistent with the grading results ( Figures 13D, E ). Both X-ray and MRI analyses confirmed that the knockdown of IRF7 alleviated IDD. Histological analysis further validated the therapeutic effects of IRF7 knockdown. It was found that NPs in the Si-NC group gradually shrank with blurred borders and were replaced by AF contents. In contrast, the NPs in the Si-IRF7 group showed significantly restored morphology ( Figure 13F ). The expression of COL2A1, MMP13, NLRP3 and IL-1β in NP tissues was accessed with immunohistochemical staining. The knockdown of IRF7 significantly restored the expression of COL2A1, while the expression of MMP13, NLRP3 and IL-1β was significantly reduced ( Figure 14 ). These findings collectively demonstrate that the knockdown of IRF7 can inhibit inflammation and reverse the metabolic imbalance of the ECM, thereby alleviating the progression of IDD.

To identify potential therapeutic drugs for IDD, 10 small molecule compounds most likely to serve as potential drugs were predicted based on the CMAP database ( Supplementary Table 15 ). Homoharringtonine, which exhibits intense toxicity, was excluded from the analysis as a potential drug. Additionally, the structure of QW-BI-011 could not be found. Therefore, eight small-molecule compounds were identified as potential drugs: Tyrphostin-AG-835, Periplocymarin, Mirin, MLN-4924, AZ-10417808, Kinetin-riboside, Securinine, and Quinoclamine. Subsequently, molecular docking analysis of potential small-molecule drugs with IRF7 was performed. The molecular docking demonstrated that the binding energies of IRF7 with these potential small molecule drugs were less than -3 kcal/mol, indicating strong binding affinity. This suggests that these small-molecule drugs may exert therapeutic effects by binding to IRF7 for IDD treatment ( Figure 15 ).