Seven expression matrix files of SLE PBMC samples were downloaded from the GEO database, namely GSE65391, GSE121239, GSE61635, GSE4588, GSE50772, GSE99967, and GSE24706. The GSE65391 dataset contained 924 SLE and 72 control samples, the GSE121239 dataset contained 292 SLE and 20 control samples, GSE61635 dataset contained 99 SLE and 30 control samples, the GSE4588 dataset contained 15 SLE and 19 control samples, the GSE50772 dataset contained 61 SLE and 20 control samples, the GSE99967 dataset contained 42 SLE and 17 control samples, and GSE24706 dataset contained 28 SLE and 20 control samples. Then, the three files GSE65391, GSE121239, and GSE61635 were pooled into a metadata cohort for the following analysis, after the batch effects were removed using the R package “SVA”. In addition, the four data sets (GSE4588, GSE50772, GSE99967, and GSE24706) were also merged as the verification data after normalization.

The DEGs were screened out using the R package “limma” based on the metadata cohort’s data set. The heat map showing the expression of DEGs was depicted using the R package. Next, the DAVID database was used to analysis the functions of the DEGs based on gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. In addition, the R package “clusterProfiler” and “DOSE” were used for the Disease ontology (DO) enrichment analyses. The gene set enrichment analysis (GSEA) was used to distinguish the most significant functional items between SLE and controls. The gene set “c2.cp.kegg.v7.0.symbols.gmt” from the Molecular Signatures Database (MSigDB) (22) was selected as the reference gene set.

Two algorithms, namely least absolute shrinkage and selection operator (LASSO) logistic regression and support vector machine-recursive feature elimination (SVM-RFE), were used to screen the potential SLE diagnostic markers. The “glmnet” R package was used to perform the LASSO analysis to identify the optimal diagnostic markers in SLE, and the SVM-RFE was used to find the optimal variables. The candidate diagnostic markers were picked by intersecting the markers identified by these two algorithms. Then, the expression of the candidate diagnostic biomarkers were verified based on the merged dataset containing GSE4588, GSE50772, GSE99967, and GSE24706.

To examine the diagnostic efficacy of the candidate markers, the receiver operating characteristic (ROC) curves were drawn based on the discovery metadata and the validation data set.

A total of 23 whole blood samples (including 13 SLE samples and 10 control samples) were collected from Shenzhen People’s Hospital. Patients diagnosed as SLE, with a SLE disease activity index (SLEDAI) score more than 5 were included. The clinical manifestations of all patients were shown in Table 1. All participants volunteered to participate in this research. This work was approved by the Ethics Committee of Shenzhen People’s Hospital (LL-KY-2021393). After collecting all of the whole blood samples, the PBMCs were isolated, and dissolved in Trizol (Beyotime, R-0016), and then stored in -80 °C.

The total RNA was isolated from the PBMC according to the manufacturer’s instructions. The transScript all-in-one first-strand cDNA synthesis superMix for qPCR (One-step gDNA removal) kit (TransGen Biotech, AT341-01) was used for the reverse transcription of mRNA. After that, the qPCR assays were conducted out with the PerfectStart Green qPCR SuperMix kit (TransGen Biotech, AQ601-02), with the following primers: GAPDH (Forward: 5-TGACTTCAACAGCGACACCCA-3. Reverse: 5-CACCCTGTTGCTGTAGCCAAA-3). ABCB1 (Forward: 5-GTCTGGACAAGCACTGAAAGATAAGA-3. Reverse: 5- CAACGGTTCGGAAGTTTTCTATTGC-3). EIF2AK2 (Forward: 5- AGAAGGCGGAGCGTGAAGTAAAAG-3. Reverse: 5- ATCCATCCCAACAGCCATTGTAGTG-3). HERC6 (Forward: 5- CCTGCCAAGCCTAAACCTGAGAAG-3. Reverse: 5- TACAGAGCCAGTGGGAAAGGAAGG-3). ID3 (Forward: 5-GGACGACATGAACCACTGCT-3. Reverse: 5-TCAGTGGCAAAAGCTCCTTT-3). IFI27 (Forward: 5-TCTGCAGTCACTGGGAGCAA-3. Reverse: 5-CCCAATGGAGCCCAGGAT-3), PLSCR1 (Forward: 5-CCTCAGTATCCACCGACAGCATTC-3. Reverse: 5-ACTGGCTGATTTGGGACAGGAAAG-3). All the primers were synthesized by the Sangon Biotech Company. Among them, GAPDH was an internal reference gene. The gene relative expressions were calculated by the 2-ΔΔCT method. A p value < 0.05 was considered statistically significant.

The CIBERSORT algorithm (23) was used to calculate the ratios of various immune cells in PBMC of SLE patients and healthy people based on the expression matrix, and the R package “vioplot” was used to visualize the proportions of 22 immune cells in SLE and control groups. The “corrplot” package was used to create the heat map displaying the quantitative correlation between distinct immune cells. In addition, the “ggplot2” R package was utilized to investigate the association between the expression of the diagnostic markers and the ratios of immune cells.

After batch effects were removed using the R package “SVA,” we combined three expression array data sets collected from the GEO database (GSE65391, GSE121239, and GSE61635) into one discovery data set, which included data from 1315 SLE patients and 122 healthy people. We used principal component analysis (PCA) to cluster all of the samples and discovered that each sample’s point was randomly distributed, indicating that the normalization was done correctly (Supplementary Figure S1). As a consequence, 284 genes were discovered to be differentially expressed in the PBMC of SLE patients versus healthy people (p value < 0.0001, fold change > 1.5) (Figure 1A). We used the heat map showing the expression of the DEGs in each healthy people and SLE patient (Figure 1B). Through an enrichment analysis of the DEGs, we found that these DEGs were enriched in some typical autoimmune-disease-relevant pathways, such as type І interferon signaling pathway, innate immune response, inflammatory response, and systemic lupus erythematosus (Figures 1C–F). In addition, we conducted the DO enrichment analyses of the DEGs, and discovered that the DEGs were primarily involved in several immune related diseases, namely bacterial infectious disease, hematopoietic system disease, human immunodeficiency virus infectious disease, and rheumatic disease, etc. (Figure 1G). Furthermore, the GSEA results showed that the DEGs in SLE patients were mainly involved with chemokine signaling pathway, complement and coagulation cascades, cytoplasmic DNA sensing pathway, RIG-I-like receptor signaling pathway, and Toll-like receptor signaling pathway (Figures 1H, I).

Then, we adopted two machine learning algorithms to screen the biomarkers for SLE, namely LASSO regression and SVM-RFE algorithm. As a result, the LASSO regression algorithm revealed 70 probable biomarkers, while the SVM-RFE approach identified 37 (Figures 2A, B). After that, we made an intersection of the 70 and 37 probable biomarkers to arrive at 14 common biomarkers (Figure 2C; Supplementary Table S1). To further confirm the reliability of these biomarkers, we verified their expressions based on the validation data set (n_SLE = 146, n_normal = 76). The result showed that 10 genes harbored the similar expression trend with statistical significance in both discovery and validation data sets, including ABCB1, EIF2AK2, HERC6, PLSCR1, ID3, IFI27, SCRN1, CD160, HSP90AB1, and PCYOX1. Among these 10 genes, we selected ABCB1, EIF2AK2, HERC6, PLSCR1, ID3, and IFI27 for the following study because the statistical differences of their expression was most significant in both the discovery and validation sets (Figures 2D–I).

Subsequently, we plotted the ROC curves for the six candidate biomarkers, and found that ABCB1, EIF2AK2, HERC6, ID3, IFI27, and PLSCR1 all had a good diagnostic effect in the discovery data set, with an AUC value 0.839, 0.912, 0.894, 0.891, 0.902, and 0.907, respectively. When we used the six markers as one single signature, the AUC value reached 0.96 (Figure 3A). On the other hand, we also verified the diagnostic effect of the six genes in the validation data set, and discovered that the AUC values of the six biomarkers, including ABCB1, EIF2AK2, HERC6, ID3, IFI27, and PLSCR1, were 0.754, 0.854, 0.81, 0.731, 0.875, and 0.851, respectively. When the six markers were combined as one single signature, the AUC value reached 0.913 (Figure 3B). The result suggests that the six genes have a good diagnostic power for SLE, and the power is higher when they are used together.

As mentioned above, SLE is highly heterogeneous among different races and regions (24). To ensure that the biomarkers we screened could be applied in the Chinese population, we further performed the RT-qPCR assays on the expression of these six markers in the PBMC of 13 Chinese SLE patients and 10 Chinese healthy controls (Supplementary Table S1). Consistently, the mRNA expression level of ABCB1 in PBMC of SLE patients was significantly lower than that of healthy controls (Figure 4A). SLE patients had considerably greater levels of IFI27 and PLSCR1 expression than healthy controls (Figures 4B, C). Meanwhile, we found that the expression differences of EIF2AK2, HERC6 and ID3 in SLE patients and healthy controls followed the same pattern as the discovery and validation cohorts (though the difference was not statistically significant) (Figures 4D–F). Thus, these results indicate that our findings are reproducible, and ABCB1, IFI27 and PLSCR1 are more likely to be SLE biomarkers in the Chinese population. Notably, ABCB1 is a novel found biomarker for SLE that has not yet to be published, to the best of our knowledge.

The onset of SLE causes changes in the proportion and function of a series of immune cells. To find the biomarkers whose expression were correlated with the proportions of immune cells, we first analyzed the ratio changes of 22 immune cells in 1315 SLE patients and 122 healthy people using the CIBERSORT algorithm. The results showed that compared with the control group, the proportions of naive B cells, CD8 T cells, CD4 memory resting T cells, CD4 memory activated T cells, and resting natural killer (NK) cells in SLE were significantly lower, while monocytes, macrophages M0, activated dendritic cells, neutrophils were higher (Figure 5A). Consistently, a number of previous studies have also demonstrated that the ratio of resting NK cells was decreased in SLE patients and the ratios of monocytes, neutrophils were increased in SLE patients (25–27).

Following that, we investigated the correlation between the ratios of the 22 types of immune cells in SLE patients, and discovered that the degrees of memory B cells and activated dendritic cells, the levels of regulatory T cells (Tregs) and memory B cells, the levels of Tregs and activated dendritic cells all had a strong positive link, respectively. Furthermore, the ratio of CD8 T cells was adversely linked with that of neutrophils, memory B cells, and monocytes (Figure 5B).

Finally, we looked at the relationship between the immune cell ratios and the expression of ABCB1, IFI27 and PLSCR1 in SLE patients. As a result, the ABCB1 expression was positively correlated with the levels of CD8 T cells, resting NK cells, CD4 memory resting T cells, naive B cells, CD4 memory activated T cells, and macrophages M2, and negatively correlated with the ratios of activated mast cells, resting mast cells, CD4 naive T cells, Tregs, memory B cells, plasma cells, macrophages M0, activated dendritic cells, neutrophils, and monocytes (Figure 5C). The IFI27 expression was positively linked with the degrees of activated dendritic cells, plasma cells, resting mast cells, activated NK cells, neutrophils, macrophages M0, memory activated cells CD4, and macrophages M1, and negatively linked with CD8 T cells, resting dendritic cells, naive B cells, and CD4 memory resting T cells (Figure 5D). The PLSCR1 expression was positively correlated with the ratios of neutrophils, activated dendritic cells, resting mast cells, macrophages M0, plasma cells, monocytes, and activated NK cells, and negatively correlated with eosinophils, CD4 naive T cells, resting dendritic cells, resting NK cells, CD8 T cells, naive B cells, CD4 memory resting T cells (Figure 5E). All in all, the ABCB1 expression was most closely linked to the ratios of CD8 T cells and monocytes, the IFI27 expression to the levels of activated dendritic cells and CD4 memory resting T cells, and the PLSCR1 expression to the degrees of neutrophils and CD4 memory resting T cells.