Our study extracted data from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Normalized microarray gene expression data of GSE49454 (whole blood, Supplementary Tables S1, S2 ), GSE61635 (whole blood, Supplementary Table S3 ), GSE110169 (whole blood, Supplementary Tables S4, S5 ), and GSE72509 (whole blood, Supplementary Table S6 ) were used as training sets to screen key m6A regulators. Out of these four datasets, only diagnostic data for GSE49454 were available. GSE50772 (peripheral blood mononuclear cells, Supplementary Table S7 ), GSE81622 (peripheral blood mononuclear cells, Supplementary Table S8 ), GSE122459 (peripheral blood mononuclear cells, Supplementary Table S9 ), GSE20864 (whole blood, Supplementary Table S10 ), GSE39088 (whole blood, Supplementary Table S11 ), and GSE156751 (plasmablast B cells, naive B cells, and memory B cells, Supplementary Table S12 ) were available from the GEO database as verification sets. The ComBat function of the R software package sva was selected to eliminate batch differences (15). The basic information of the datasets selected is shown in Supplementary Table S13 .

To screen m6A regulators, all samples from 10 datasets were analyzed using the empirical Bayesian method of the limma package of R. The criterion for screening m6A regulators was an adjusted p-value <0.05. The common m6A-mediated genes overlapped according to a Venn plot. Enrichment analyses of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were conducted using Metascape (http://metascape.org/gp/index.html#/main/step1) and FunRich (functional enrichment analysis tool) (http://www.funrich.org/), respectively. The coexpression networks of the m6A regulators were constructed in GENEMANIA software (http://genemania.org/search/).

The comparative toxicogenomics database (CTD, http://ctdbase.org/) integrates information, including gene–disease relationships, to explore the underlying pathogenesis of diseases (16). The association between potential crucial m6A genes and SLE risk was analyzed using the CTD data.

The interaction network of microRNAs (miRNAs) and transcription factors (TFs) associated with the key m6A regulators was predicted by Networkanalyst (https://www.networkanalyst.ca/).

The “corrplot” R (17) package was used to perform the correlation analysis of m6A regulators or immune cell infiltration. When the correlation coefficient is >0.5, the correlation is positive. When the correlation coefficient is <0.5, the correlation is negative.

RNA-seq data of 23 m6A regulators, including 8 “writers” (METTL14, METTL3, WTAP, METTL16, VIRMA, RBM15B, ZC3H13, and RBM15), 13 “readers” (YTHDF1, HNRNPA2B1, YTHDF2, YTHDC1, IGF2BP1, RBMX, HNRNPC, LRPPRC, YTHDC2, IGF2BP2, IGF2BP3, YTHDF3, and FMR1), and 2 “erasers” (ALKBH5 and FTO), were obtained from the datasets. Least absolute shrinkage and selection operator (LASSO) regression was used to construct the m6A diagnostic signature by feature selection and dimension reduction (18). Cross-validation was used for lambda parameter tuning. The cutoff value was decided by rank statistics. Unsupervised clustering analysis was used to identify different m6A modification patterns. The cluster numbers and robustness were evaluated as described in a previous study (19). The R package “ConsensusClusterPlus” was used for classification (20). PCA was used to validate the signature score. The m6A gene expression and immune cell infiltration abundance scores among the distinct modification patterns were compared. The enrichment scores of immune cell abundance and immune reactivity in healthy and SLE samples were compared by Wilcoxon assay.

In this study, the “ConsensusClusterPlus” R (20) package was used to cluster samples based on gene expression, and 1,000 iterations were performed to ensure the stability of classification. We used the “CIBERSORT” package R (21) to quantify the levels of infiltration of different immune cells in different datasets. Different immune cell subtypes were grouped and utilized to compare the discrepancies of different clusters. Dimension reduction in ICI gene clusters was based on the Boruta algorithm and the signature score using PCA. ESTIMATE and CIBERSORT algorithms were used for the analysis of immune status, such as immune and stromal scores (21, 22). For the ICI score, the following method was used (23): ICI score = ∑PCI_A − ∑PCI_B.

We downloaded gene expression and clinical characteristic data to perform univariate and multivariate Cox analyses to identify whether the expression of m6A regulators or immune genes was related to the diagnosis of SLE patients. Several clinicopathological parameters were included. The m6A regulators or immune genes were identified by univariate logistic regression (cutoff criterion: p-value <0.05). Multivariate logistical regression was applied to verify m6A regulators related to the SLE classifier.

Gene signature univariate Cox regression analysis (p < 0.001 was considered statistically significant) and LASSO‐penalized Cox regression analysis (24) were used to verify the diagnostic gene signature. Diagnostic analysis of m6A immune-related genes was identified using the “survminer” and “survival” R packages (p < 0.05 was statistically significant) (25). In particular, the “survminer” R package was used to divide patients into high‐ and low‐risk groups according to the optimal cutoff value. The “survivalROC” (26) R package was used to determine the time‐dependent diagnostic value of the gene signature (log‐rank p < 0.05 was considered significant).

The list of infiltrating immune cell gene sets was derived from a previous study (19), and the immune gene sets were obtained from the ImmPort database (http://www.immport.org) (27). Univariate Cox regression analysis and LASSO-penalized Cox regression analysis were applied to verify the diagnostic signature in the training dataset. The risk score was calculated as follows: risk score = SExpi*bi. Expi represents each gene expression, and bi represents the coefficient of each gene.

Weighted gene coexpression network analysis (WGCNA) was performed using the WGCNA R package (28) to build a gene coexpression network, and key genes related to the significant modules were identified. Gene modules with |correlation coefficient| >0.5 were considered strongly correlated modules, and significantly diagnosis-associated hub immune genes were selected for further analysis (p < 0.05, log-rank test). Strong diagnosis-correlated gene modules were used as an input in the LASSO regression analysis.

The GO (29) R package was used to display the results of the GO and KEGG analyses. All differentially expressed m6A immune genes were selected for GO and KEGG pathway analyses. GO analysis covers three domains: biological process (BP), cellular component (CC), and molecular function (MF). The importance of the GO and KEGG pathways was identified [false discovery rate (FDR) <0.05 denotes the significance of the p value].

The Gene Set Variation Analysis (GSVA) package in R (30) was used to evaluate the most significantly enriched molecular pathways. The gene sets of “c2.cp.kegg.v7.0.symbols” and “h.all.v7.0.symbols’ were downloaded from the Molecular Signatures Database (MSigDB, version 6.0) database for running GSVA analysis. Differential analysis of the enrichment scores of KEGG pathways between SLE patients and healthy controls was performed using the R package limma (31).

The significantly enriched pathways in samples with ICI score-high or ICI score-low were identified by GSEA according to the MSigDB. The number of infiltrating immune cells and the activity of the immune response were estimated by ssGSEA, and the enrichment fraction represents the absolute enrichment degree of a gene set in each sample (32).

R software (https://www.r-project.org/, version 4.1) and corresponding R packages (http://www.bioconductor.org/) were utilized for statistical analysis. The results are presented as the mean ± SD (p < 0.05 indicated statistical significance). Student’s t-test was applied for continuous variable analysis. Spearman’s correlation was adopted to analyze correlations. Kaplan–Meier analysis was performed, and diagnostic values were compared by the log-rank test.

Based on the high-throughput data analysis, four whole-blood datasets were used as training sets to screen the key m6A regulators based on the expression levels of m6A genes. The results showed that 4 m6A regulators were significantly differentially expressed in the GSE49454 dataset ( Figure 1A ), 17 m6A regulators were significantly differentially expressed in the GSE61635 dataset ( Figure 1B ), 13 m6A regulators were significantly differentially expressed in the GSE110169 dataset ( Figure 1C ), and 12 m6A regulators were significantly differentially expressed in the GSE72509 dataset ( Figure 1D ). The principal component analysis of GSE49454, GSE61635, GSE110169 and GSE72509 datasets were performed by unsupervised clustering based on the expression of genes ( Figure 1E ). The Venn diagram showed that IGFBP3 overlapped in the four datasets ( Figure 1F ), which suggests that it plays an important role in SLE patients.

GO and KEGG pathway analyses revealed the functions of the 19 m6A regulators. The GO results showed that they were mainly enriched in the regulation of mRNA metabolic processes, mRNA processing, and related processes ( Figure 2A ). The KEGG results showed that these regulators were mainly enriched in the regulation of insulin-like growth factor (IGF) activity by the insulin-like growth factor-binding protein (IGFBP) pathway ( Figure 2B ). In addition, the gene–gene interaction network for 19 m6A regulators was constructed by the GeneMANIA database. The hub nodes represent genes that were significantly correlated with 19 m6A regulators ( Figure 2C ). CTD was employed to explore the interaction between potential crucial genes and SLE. Inference scores in CTD reflected the association between disease and 19 m6A regulators. The interaction results showed that IGFBP3 has a higher inference score with SLE ( Figure 2D ).

To verify whether IGFBP3 could act as a key target in SLE, based on high-throughput analysis, six datasets were selected as verification sets to screen significantly differentially expressed m6A regulators after the ChIP results were normalized. Among the three PBMC datasets, nine, six, and five m6A regulators were significantly differentially expressed in the GSE50772 dataset ( Supplementary Figure S1A ), the GSE81622 dataset ( Supplementary Figure S1B ), and the GSE122459 dataset ( Supplementary Figure S1C ), respectively. In addition, the m6A regulators in two whole-blood microarray datasets were screened. Four and 15 m6A regulators were significantly differentially expressed in the GSE20864 dataset ( Supplementary Figure S1D ) and the GSE39088 dataset ( Supplementary Figure S1E ), respectively. Four m6A regulators were significantly differentially expressed in the B-cell dataset of GSE156751 ( Supplementary Figure S1E ). The Venn diagram showed that IGFBP3 overlapped in the six verification datasets ( Figure 3A ). Furthermore, IGFBP3 overlapped in the three subsets [GSE156751 (plasmablast B cells), GSE156751 (memory B cells), GSE156751 (naive B cells)] of the GSE156751 dataset ( Figure 3B ). These results also indicated that IGFBP3 could act as a key target in B cells to further explore the detailed molecular mechanism of SLE. Then, GSVA was applied to explore the molecular pathways and underlying mechanisms in B-cell subsets between SLE patients and healthy controls. The top differentially enriched molecular pathways were identified ( Figures 3C–E ). The results showed that olfactory transduction and neuroactive ligand receptor interaction pathways were positively correlated with SLE. DNA replication, RNA polymerase, neurotrophin signaling, spliceosome, etc. pathways were negatively correlated with SLE.

Univariate Cox regression analysis was used to identify diagnosis-related m6A regulators. Seven genes (HNRNPA2B1, IGFBP1, WTAP, YTHDC1, YTHDC2, YTHDF1, and YTHDF3) were related to increased risk (HRs > 1), including 14 protective genes with HRs <1 ( Figure 4A ). LASSO Cox regression analysis constructed the signature according to the optimum λ value ( Figures 4B, C ). Receiver operating characteristic (ROC) analysis was used to evaluate the sensitivity and specificity for the diagnostic model between the low- and high-risk groups ( Figure 4D ). Furthermore, a heatmap of clinical features for the GSE49454 cohort is shown ( Figure 4E ), and the diagnostic status, low c3, low c4, and age of patients were diversely distributed between the low- and high-risk subgroups (p < 0.05). Meanwhile, we performed multivariate Cox analysis of IGFBP3 for some variables (age, sex, low c3, and low c4) in GSE49454. The multivariate Cox analysis indicated that IGFBP3 (HR = 1.51, p = 0.0089, 95% CI = 0.94–2.42) and age (HR = 0.96, p < 0.001, 95% CI = 0.94–0.97) were more accurate than the other factors ( Figure 4F ). Next, we obtained two m6A risk genes from the intersection of the six m6A risk model genes between the low- and high-risk subgroups and four m6A regulators in GSE49454 ( Figure 4G ). For the two common m6A risk genes (IGFBP2 and IGFBP3), a TF–miRNA gene interaction graph was built using NetworkAnalyst ( Figure 4H ). The network contains a total of 26 TF genes and 27 miRNAs. This interaction might be the reason for regulating the expression of the m6A risk model genes.

Consensus cluster analysis was used to develop the molecular classification of SLE patients according to gene expression. All SLE samples were initially divided into different k (k = 2–9) groups. The cumulative distribution function (CDF) curves showed that the optimal number of clusters was k = 3. Consensus clustering was independently conducted in the GSE49454 ( Figures 5A, B ) and GSE110169 ( Figures 5C, D ) datasets, as shown in the heatmap, in which SLE patients were categorized into three m6A and gene clusters. The two-dimensional principal component plot showed that the clusters of samples were separated from each other ( Figures 5A, C ).

Infiltrating immune cells were detected to understand the differences in the characteristics of immune infiltration among m6A clusters. Immunocytes differed in GSE49454 and GSE110169. Activated.CD8.T.cellna, Activated.dendritic.cellna, macrophagena, type.1.T.helper.cellna, T.follicular.helper.cellna, mast.cellna, monocytena, neutrophilna, and type.17.T.helper.cellna were significantly different in GSE49454 ( Figure 6A ), while Activated.CD4.T.cellna, Activated.CD8.T.cellna, Activated.dendritic.cellna, immature.B.cellna, immature.dendritic.cellna, CD56bright.natural.killer.cellna, mast.cellna, natural.killer.cellna, plasmacytoid.dendritic.cellna, type.1.T.helper.cellna, and type.2.T.helper.cellna were markedly different in GSE110169 ( Figure 6B ). Furthermore, the expression differences of m6A regulators between three gene clusters in GSE49454 and GSE110169 were evaluated. Sixteen m6A regulators were significantly different in GSE49454 ( Figure 6C ), and 21 m6A regulators were significantly different in GSE110169 ( Figure 6D ). Moreover, there were obvious differences in m6A score in the molecular subtypes (gene clusters and m6A clusters) of SLE in the GSE49454 ( Figures 6E, F ) and GSE110169 ( Figures 6G, H ) datasets.

Unsupervised consensus clustering was used to explore the molecular classification of SLE patients based on the expression patterns of immune cell infiltration (ICI). According to the relative change in the area under the CDF curve, the optimal number of clusters was determined to be three (k = 3). Hence, all SLE patients were categorized into three groups, which were termed ICI clusters 1–3. The differences in immune infiltration characteristics among the three m6A clusters and two gene clusters are shown by heatmaps and boxplots in GSE110169 ( Figures 7A–D ). The results revealed that naive CD4⁺ T cells, CD8⁺ T cells, activated memory CD4⁺ T cells, monocytes, gamma delta T cells, M2 macrophages, neutrophils, and immune scores were significantly different in the ICI clusters of GSE110169 ( Figures 7A, B ). Gamma delta T cells and monocytes were significantly different in the gene clusters of GSE110169 ( Figures 7C, D ). Consistently, to explore the differences in immune infiltration characteristics among three ICI clusters and two gene clusters, infiltrating immunocytes were also evaluated. Many immunocytes differed between the two patterns in GSE49454 ( Figures 8A–D ). Memory B cells, naive CD4⁺ T cells, CD8⁺ T cells, activated memory CD4⁺ T cells, monocytes, M0 macrophages, activated dendritic cells, neutrophils, stromal scores, and neutrophils were significantly different in the ICI clusters of GSE49454 ( Figures 8A, B ). Memory B cells, CD8⁺ T cells, activated dendritic cells, stromal scores, and neutrophils were significantly different in the gene clusters of GSE49454 ( Figures 8C, D ).

A complete immune gene map was constructed. The related gene modules of different traits were identified by WGCNA. Two gene modules were identified to match different clinical characteristics, with the diagnosis closely associated with the genes in the blue module ( Figure 9A ). Univariate logistic regression was applied to identify immune genes associated with SLE, and five immune regulators (TAP2, DDX58, IDO1, CD14, and FGFRL1) were related to SLE ( Figure 9B ). Three key immune genes (FGFRL1, IDO1, and CD14) were selected by LASSO regression ( Figure 9B ). GO and KEGG enrichment analyses were used to investigate the biological pathways of the immune genes. The results revealed that they are mainly involved in neutrophil activation ( Figure 9C ). The biological processes involved in immune genes are closely related to T-cell activation pathways ( Figure 9D ).

To screen for key genes, a total of 6,148 common genes were obtained from the intersection of clustered genes in GSE110169 ( Figure 10A ), and GO and KEGG enrichment analyses demonstrated that they mainly participated in the process of detection of chemical stimuli involved in sensory perception and olfactory transduction pathways ( Figures 10B, C ). GO and KEGG pathway analyses were also used to evaluate the key biological pathways of immune genes in GSE110169. GO enrichment analysis showed that they are mainly involved in the processes of neutrophil degranulation and neutrophil activation during the immune response ( Figure 10D ). The biological pathways in which the immune genes take part were markedly related to CENP-A-containing nucleosome assembly, DNA replication-independent nucleosome assembly, etc. ( Figure 10E ).

We used univariate Cox proportional hazards model analysis to build a risk score model based on the median-risk score, which divided SLE patients and healthy samples into two groups: high risk and low risk in the GSE49454 dataset ( Figures 11A–C ). By LASSO regression analysis, 13 differentially expressed immune genes (PDIA3, RAET1L, LCN8, DEFB107B, ISG20, IDO1, SOCS1, CD14, IL1RN, FGFRL1, GCGR, NGFR, and CASP3) were screened out. Furthermore, we quantified the enrichment levels of the 13 genes in the two groups, and the results are shown in a heatmap ( Figures 11A–C ) and a boxplot ( Figures 11D–F ). The average AUC values for the sensitivity and specificity of the risk score and clinical factors in the three sets reached 0.733, 0.725, and 0.741. ( Figures 11D–F ). Kaplan–Meier curve analyses ( Figures 11G–I ) showed the diagnostic probability between high- and low-risk patients in the GSE49454 cohort (p < 0.001).

The correlation coefficient heatmap provides an intuitive understanding of the state of immune cell interaction of GSE110169 ( Figure 12A ). CIBERSORT was used to calculate the immune cell expression of SLE subtypes in the GSE49454 dataset. The Wilcoxon test was used to compare the distribution of immune cells between high- and low-risk patients. We found that the differences in most immune cells were insignificant in the high-risk group compared with the low-risk group ( Figures 12B–D ). Only naive CD4⁺ T cells (p = 0.039) and M0 macrophages (p = 0.006) were higher in the low-risk group. In contrast, resting dendritic cells (p < 0.001) and resting memory CD4⁺ T cells (p = 0.025) were higher in the high-risk group ( Figure 12D ).

The samples were unsupervised clustered using R packages. The graphs of unsupervised clustering of GSE49454 are shown in Figures 5 , 8 . Significant differences in dependent clusters on diagnosis were revealed (log-rank test, p < 0.001) ( Figures 13A–C ). To examine the effectiveness of the m6A score in predicting clinical features, subjects were divided into high or low m6A score groups. The results showed that the number of SLE patients with high m6A scores was higher than that of patients with low m6A scores ( Figure 13D ). The comparison showed that the m6A scores of the SLE patients were higher than those of the healthy samples ( Figure 13E ). The Wilcoxon test showed that IGFBP3 was significantly overexpressed in the high ICI group ( Figure 13F ). The correlation coefficient heatmap revealed the immune cell interaction ( Figure 13G ). By GSEA, the remarkable enrichment of KEGG pathways delineated biological pathways or processes correlated with the ICI score in the GSE49454 cohort. Five representative KEGG pathways with high and low ICI scores were identified. The results are as follows: glycine serine and threonine metabolism, renin–angiotensin system, tight junction, non-small cell lung cancer, and melanoma signaling pathways were associated with a high ICI score, while enrichment plots demonstrated that basal transcription factors, lipid metabolism, proteasome, RNA polymerase, and ribosome pathway were associated with a low ICI score ( Figure 13H ). The relationships of gene cluster, ICI score, diagnosis, and the relationship of m6A cluster, gene cluster, m6A score, and diagnosis were also shown ( Figures 13I, J ).

The AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk/) was used to predict the protein structures of IGFBP3, FGFRL1, CD14, and IDO1 ( Figures 14A–D ). We also compared the representative CD14 and IDO1 structures in the Protein Data Bank (PDB) database, which contains experimental structural data. IGFBP3 and FGFFRL1 were not included in PDB. We screened the drug-target complex from DrugBank (https://go.drugbank.com/) according to each target’s directly interacting ligands.