We retrieved datasets related to SD and their corresponding control groups from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) (12). Human RNA expression data from 17 individuals with SD and 25 healthy controls were obtained from the GSE208668 dataset, derived from the GPL10904 platform. For the validation of the subsequently identified hub genes, we utilized the GSE240851, GSE56931, and GSE98582 datasets, which were sourced from the GPL24676, GPL10379, and GPL6244 platforms, respectively. These datasets included varying numbers of SD patient and control samples. Additionally, the GSE165041 dataset, generated on the GPL18573 platform, comprised microRNA expression data from 10 SD patients and an equal number of healthy controls. All datasets were normalized using the “limma” package within the R software environment, version 4.1.2 (13).

We utilized the limma package in R software to process the normalized datasets GSE240851 and GSE208668, aiming to identify DEGs. To ascertain statistical significance, DEGs were determined based on |log2fold change| ≥ 0.58 and a false discovery rate <0.05. Subsequently, we employed the ggplot2 package to create volcano plots for the DEGs. Moreover, we selected the top 20 genes, ranked by |log2fold change|, to build heatmaps as part of our analysis.

The “c2.all.v2023.1.Hs.entrez” and “c5.all.v2023.1.Hs.entrez” datasets, which serve as reference gene sets, were downloaded from the GSEA official website originating from the MSigDB database (14). Gene set enrichment analysis (GSEA) was then performed using the “clusterProfiler” package in R software (15). Subsequently, the analysis results were visualized utilizing the “enrichplot” package in R software.

We employed the ImmuCellAI and ssGSEA methodologies to estimate the abundance of 24 and 28 distinct types of immune cells, respectively, in the tissues of SD patients and a normal population, thereby enabling precise delineation of immune cell profile differences between the two groups (16, 17). Furthermore, we conducted Spearman correlation analysis to elucidate the interrelationships among the distributions of various immune cells.

We downloaded the H: hallmark gene sets from the GSEA official website to probe the functional disparities between SD patients and the normal cohort. With these gene sets, we applied the ssGSEA algorithm to evaluate 50 gene sets, aiming to discern potential variances between the two groups. Subsequently, we utilized the Mantel algorithm to analyze the correlations among these gene sets (18). This approach allowed us to investigate the differences in gene expression profiles between the two cohorts and to identify potential molecular pathways associated with the disease.

For protein interaction analysis, we utilized the STRING 4 online platform and specifically selected PPI pairs with a confidence score greater than 0.40. Subsequently, we employed the Cytoscape V3.9.0 software for visualizing the PPI network (19). Within the network, the significance of each node was determined by calculating their Degree values using the CytoHubba plugin. This analysis allowed us to identify the top 20 pivotal genes based on their ranking of importance (20).

For gene selection, we utilized the Random Forest (RF) algorithm, a binary tree-based recursive partitioning method. The “randomForest” package in R was used, with parameters set to ntree=1000, mtry=3, and importance=true (21). Employing the Gini index as the primary assessment criterion, the Random Forest algorithm was used to rank the DEGs, and the top 20 genes with a significance value greater than 3 were earmarked for further analysis.

We utilized the SVM-RFE (Recursive Feature Elimination) approach to optimize the predictive model by minimizing the number of feature vectors produced by the SVM. This approach, being an effective binary classification tool, operates by constructing a classification hyperplane to delineate decision boundaries. In order to enhance the algorithm’s precision, we configured parameters to method=repeatedcv and repeats=10 in the R package “1071”, and employed ten-fold cross-validation. This approach aimed at augmenting the algorithm’s precision (22).

To develop the LASSO model (23), we integrated genes identified through CytoHubba, RF, and SVM algorithms. This method effectively enabled the identification of critical hub genes for diagnosing SD. We then employed the Logistic regression approach to investigate pivotal factors associated with SD, ultimately constructing a simplified model. Subsequently, we evaluated the classification performance of the model using the Receiver Operating Characteristic (ROC) curve and the corresponding Area Under the Curve (AUC).

In order to evaluate the robustness and general applicability of our developed diagnostic model, we computed the AUC of the ROC curve for the model using three distinct datasets: GSE240851, GSE56931, and GSE98582. This procedure aimed to ascertain the model’s performance across diverse datasets. It was essential to assess its efficacy in different contexts and ensure its capability to perform consistently across varied data sources.

We utilized the Limma package in R to conduct a differential analysis of the expression matrix from GSE165041, aiming to identify miRNAs that were differentially expressed (DEmiRNAs). Our criteria for significance were miRNAs exhibiting a false discovery rate <0.05 and |log2fold change| >0.5. Moreover, we employed the miRNet database to investigate potential miRNAs associated with differentially expressed genes (DEGs), in order to gain further insight into their involvement in SD (24).

We obtained data on drugs related to core genes from the Drug-Gene Interaction database (DGIdb) (25). We used the “ggplot2” package in R to create bar charts showing interaction scores, visually indicating the intensity of interactions between various drugs and core genes.

The SD model, following the methodology outlined by Alkadhi and Alhaider (26), was established in 10-month-old male mice. The experimental mice were subjected to SD treatment for a duration of 8 weeks by being placed on a small fixed platform encircled by water, with access only to water and food. Meanwhile, for the control group, another set of mice was housed in a comfortable environment with a 12-hour light/dark cycle and unrestricted access to water and food. All mice received standard pellet feed, with the daily quantity of feed being maintained uniformly across all groups. This experimental design aimed to simulate the effects of SD on physiological functions, thus laying the groundwork for subsequent investigations into gene expression and drug treatment efficacy.

The cortex in SD mice was used for total RNA extraction, employing the TransZol Up Plus RNA Kit (TransGEN, Beijing, China) (27). The RNA concentration and quality were then assessed using the Nanodrop Spectrophotometer (Termo Scientifc, Waltham, MA, USA). Subsequently, reverse transcription was carried out using the TransScript^® One-Step gDNA Removal and cDNA Synthesis SuperMix (AT311, TransGEN, Beijing, China). Amplification was monitored with the ChamQ Universal SYBR qPCR Master Mix (Novozymes Q711) and a QuantStudio™ 5 Real-Time PCR System (Thermo Fisher Scientific). The internal reference was β-actin, and the relative gene expression was determined using the 2-ΔΔCT formula. Detailed primer sequences can be found in Table 1 .

Western blotting was performed for IPO9, RAP2A, DDX17, and GAPDH, as described previously ( Table 2 ). Enhanced chemiluminescence reagents were used to detect protein expression, and quantitative analysis was conducted using Image J software (28).

The normality of the data was assessed using the Shapiro–Wilk test. The t-test was utilized to compare the data between two groups, while comparisons among multiple groups were conducted using one-way ANOVA, followed by either the LSD post hoc test or Tukey’s post hoc test. The statistical analyses were conducted using R software version 4.2.1 and SPSS 25, with significance defined as P < 0.05. Furthermore, all experiments were independently repeated at least three times to ensure the validity of the results.

In the GSE208668 and GSE240851 datasets, a differential expression analysis was conducted on SD samples and normal control samples. The analysis identified 5964 DEGs in the GSE208668 dataset and 1375 DEGs in the GSE240851 dataset. Utilizing these DEGs, volcano plots ( Figures 1A, B ) were constructed, and the top 20 genes with the highest |log2fold change| values were selected for heatmap generation ( Figures 1C, D ). These results demonstrate the effectiveness of DEGs in distinguishing between the SD group and the normal control group.

The differential gene analysis of GSE208668, using MSigDB’s C2 and C5 gene sets, unveiled substantial differential expression in gene clusters relevant to immune cell functions, stress responses, and nervous system activities. Specifically, within the C2 gene set, activation of the gene cluster associated with nervous system development was observed, while the gene cluster linked to the B lymphocyte network was suppressed. In the C5 gene set, gene clusters involved in stress response regulation exhibited a downregulated trend, indicating a potential imbalance in environmental stress response regulation. Additionally, gene clusters associated with spinal cord development demonstrated differential expression. Moreover, immune-related gene clusters, including those implicated in autoimmunity and T-cell deficiency, exhibited alterations in expression levels ( Figures 1E, F ). Further comprehensive details can be found in Supplementary 1 .

The results of the GSEA analysis demonstrated significant enrichment of terms in KEGG pathways closely associated with neurodegenerative diseases and metabolic pathways, such as Alzheimer’s disease, leishmaniasis infection, lysosomal function, oxidative phosphorylation, Parkinson’s disease, and viral myocarditis. Notably, major signaling factors involved in immune regulation, including IL-4, IL-8, and IL-12, were also found to be significantly enriched in various immune system pathways. Furthermore, signaling pathways related to neural signaling, including neurotrophic factor signaling, B cell receptor signaling, and tumor necrosis factor α signaling, demonstrated significant enrichment Supplementary 2 .

After integrating the GSE208668 and GSE240851 datasets and mitigating batch effects, we observed significant discrepancies in the distribution of various immune cell types between patients with systemic lupus erythematosus (SD) and the normal control group. This was confirmed by employing two distinct immune infiltration scoring methods ( Figures 2A–D ). Notably, the infiltration scores for Activated CD4+ T cells and Activated CD8+ T cells were substantially higher in patients with SD compared to the normal group. Moreover, central memory CD4 T cells, central memory CD8 T cells, and natural killer T cells exhibited an increasing trend in SD patients. In contrast, the infiltration scores of Activated B cells, Immature B cells, Th17 cells, and Th2 cells demonstrated a declining trend in patients with SD relative to the normal group. Correlation analysis based on both algorithms indicated a level of concordance. Specifically, CD8+ T cells and memory T cells consistently showed a strong positive correlation, while neutrophils and regulatory T cells exhibited a negative or non-significant correlation ( Figure 3A ). For a detailed correlation analysis of the ImmuCellAI algorithm scores, please refer to Supplementary 3 .

The heatmap in Figure 3B displays the scores for 50 gene sets, from which we chose 12 related to signaling pathways and 22 associated with physiological functions for further investigation. A comparative analysis of the gene set scores between the groups revealed that in the SD group, median scores for gene sets such as androgen response, apoptosis, complement, hypoxia, and inflammatory response were significantly higher, while the median score for the myogenesis gene set was notably lower compared to the normal group. Additionally, significant disparities were observed in signaling pathways such as il2 stat5 signaling, il6 jak stat3 signaling, the p53 pathway, and tgf beta signaling ( Figures 3C, D ).

We then focused on 12 gene sets with significant differences in physiological functions between the groups for further correlation analysis. These analyses revealed that, except for the inflammatory response gene set, other gene sets showed significant correlations with immune infiltration and signaling pathway scores. Moreover, the scores within these 12 biological function gene sets also exhibited a high degree of correlation ( Supplementary 4 ). For detailed results on pathway scoring and physiological function scoring correlations, please refer to Supplementary 5 .

We first identified 138 co-expressed differential genes by intersecting the DEGs from the GSE208668 and GSE240851 datasets ( Figure 4A ). Subsequently, 52 genes were selected for network visualization in Cytoscape using the STRING database for PPI analysis with a medium confidence threshold of 0.4 ( Figure 4B ). To identify the top 20 key genes within the PPI network, the CytoHubba plugin in Cytoscape was employed ( Figure 4C ).

The assessment of hub genes started with the use of the GSE208668 dataset and the RF algorithm, which led to the identification of 138 co-expressed differential genes. The model’s performance was then evaluated across a spectrum of tree quantities, showing a gradual decrease in the error rate as the number of trees increased, ultimately indicating enhanced stability and reliability ( Figure 5A ). Following this, the feature importance was assessed through the calculation of the percentage decrease in mean impurity (%IncMSE), which resulted in the identification of the top 20 genes that had the most significant contributions to the model’s prediction ( Figure 5B ). Subsequently, the SVM algorithm was utilized to determine the hierarchical importance of the 138 genes ( Figure 5C ). The model’s predictive efficacy was then measured using the Root Mean Square Error (RMSE) through cross-validation methodologies, as shown in Figure 5D . The results indicated that the model demonstrated its strongest predictive capability when N was set to 22, leading to the selection of these 22 genes for further analysis.

After identifying 40 key genes through the amalgamation of genes identified via CytoHubba, Random Forest, and SVM ( Supplementary 6 for details), LASSO regression was utilized to select diagnostic genes. This process involved adjusting the regularization parameter λ and observing its effects on coefficient estimation ( Figure 5E ). The optimal λ value was determined through the ten-fold cross-validation method ( Figure 5F ). In the final analysis, diagnostic model construction involved the selection of IPO9, RAP2A, DDX17, MBNL2, PIK3AP1, and ZNF385A.

When compared to the normal group, these 6 hub genes demonstrated significant disparities ( Figure 6A ) and significant correlations with scores of biological function genes ( Figure 6B ). The interactions between hub genes and other functional scores are presented in Supplementary 7 . A nomogram was constructed based on these 6 genes, with each gene correlating with a specific scoring criterion ( Figure 6C ). The calibration curve of the nomogram indicated commendable predictive performance of the model ( Figure 6D ). ROC curve analysis further demonstrated the substantial diagnostic value of these genes with an overall AUC value of 0.916 ( Figure 6E ). To corroborate the accuracy of the model, three independent datasets were utilized in logistic regression analysis. The AUC values for GSE240851, GSE98582, and GSE56931 were 0.835, 0.863, and 0.751, respectively, demonstrating the stability and reliability of the model in diagnosing SD ( Figures 6F–H ).

The analysis of the GSE165041 dataset revealed 9 upregulated DEmiRNAs that distinguished the SD group from the normal group, as evidenced in Figure 7A . Subsequently, an exploration of the association between these DEmiRNAs and the previously identified 138 DEGs was conducted using the miRnet database, resulting in the identification of 5 miRNAs with potential as therapeutic targets, as depicted in Figure 7B . Following this, a validation of the potential therapeutic efficacy of the 5 identified miRNAs and their corresponding key mRNAs was performed, as shown in Figure 7C . Subsequently, compounds with high interaction scores with the key genes were identified utilizing the DGIdb, as illustrated in Figure 7D . These findings indicate that the identified miRNAs and compounds might serve as potential therapeutic agents for the treatment of SD. Nevertheless, it is essential to note that further research is necessary to substantiate their efficacy and safety.

RT-qPCR was used to scrutinize the expression levels of 6 hub genes in the tissue. The mRNA expression levels of 6 genes in the SD group showed a significant elevation relative to the normal group, as illustrated in Figure 8A . Moreover, Western blotting analysis revealed significant elevations in the protein levels of genes such as IPO9, RAP2A, and DDX17 within the cortical regions of mice suffering from SD, as depicted in Figures 8B, C . These findings align with previous bioinformatics analysis outcomes, providing further support for the correlation of expression levels among these 6 hub genes.