“Systemic Lupus Erythematosus” and “RNA-seq” were used as the keywords for searching the GEO database. The gene expression datasets of PBMC from freshly isolated healthy controls and SLE patients were downloaded from the GEO database (GSE162828 and GSE169080), the platforms used were GPL24676, and GPL20795. GSE162828 included 10 samples of peripheral blood mononuclear cells and was divided into the SLE group (5 samples) and healthy controls group (5 samples). GSE169080 included seven samples of peripheral blood mononuclear cells and was divided into SLE group (4 samples) and healthy controls group (3 samples) (Clough and Barrett, 2016; Cheng et al., 2021; Jiang et al., 2021).

R package DESeq2 (1.26.0) was used for the analysis of the original datasets (Love et al., 2014). |log FC| > 1 and p. adj <0.05 were defined as the cutoff values for further analysis of DEGs. Volcano and heatmap were constructed by R package ggplot2. Venn plot (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to draw the intersection of two databases.

Protein-protein interaction (PPI) network analysis helps to study the molecular mechanism of diseases from a systematic perspective and discover new drug targets (Wu et al., 2019). STRING (https://string-db.org/) is a database covering more than 5,000 organisms with known and predicted protein-protein interactions, providing direct (physical) and indirect (functional) associations (Szklarczyk et al., 2017). We used String (https://string-db.org/) to generate biological networks for proteins, and the results were analyzed by Cytoscape (Shannon et al., 2003; Szklarczyk et al., 2017).

Gene Ontology (GO) is an ontology widely used in the field of bioinformatics, which covers three aspects of biology: biological process (BP), cellular component (CC), and molecular function (MF) (Thomas, 2017). Kyoto Encyclopedia of Genes and Genomes (KEGG) is a biological system advanced function and utility database based on molecular-level information from genome sequencing and other high-throughput experimental technologies (Kanehisa et al., 2017). In this study, R package clusterProfiler was used to perform GO functional annotation and KEGG pathway enrichment analysis for DEGs (Yu et al., 2012).

To obtain abnormally expressed genes in SLE patients, we separately analyzed the differential expression of two GEO datasets (GSE162828 and GSE169080). As shown in Figure 1A, there were 4,713 DEGs, including 2,717 up-regulated and 1,996 down-regulated in the GSE162828 dataset. In the GSE169080 dataset, there were 2,473 DEGs, including 1,552 up-regulated and 921 down-regulated (Figure 1B). In both datasets, the number of up-regulated DEGs was more than the number of down-regulated DEGs (Figure 1C). In the GSE162828 dataset, the up-regulated DEGs accounted for 56.7% of all DEGs. At the same time, the up-regulated DEGs accounted for 62.8% of all DEGs in the GSE169080 dataset. The trends in the two datasets were roughly the same.

In addition, the heatmap showed that DEGs can group samples by sample type, namely SLE patients (SLE) and healthy controls (Norm) (Figures 1D,E). These genes were highly concordant within groups. The expression level of these genes between SLE patients and healthy controls exhibited a large difference in both databases.

Due to the heterogeneity between different datasets, the analysis results of different datasets may have certain differences (Ying et al., 2020). The gene expression in different samples may be different (Bao et al., 2021). To avoid this problem, integrating multiple datasets and a large number of samples help obtain more solid results (Kou et al., 2020). In this study, we integrated DEGs from two datasets to obtain common DEGs.

The Venn diagram showed that 790 DEGs were shared between the two datasets (Figure 2A). They accounted for 16.8% (GSE162828) and 31.9% (GSE169080) of the two datasets, respectively. There were 3,923 DEGs only in the GSE162828 dataset, and 1,683 DEGs only in GSE169080 dataset. This may be caused by different sequencing technologies and sample heterogeneity.

We defined these 790 DEGs as common DEGs. To further explore the distribution of common DEGs on the chromosomes, we had made statistics on the chromosomal locations of these genes. As shown in Figure 2B, we found that these genes were distributed on every chromosome. Most of these genes were distributed on chromosome 19. On the contrary, they were only 6 DEGs on chromosome 13.

Proteins usually perform biological functions in concert. It has been shown that there is a close relationship between Protein-Protein Interaction (PPI) and the biological functions of gene/protein clusters (Li H. et al., 2019). To further analyze the correlation between common DEGs, STRING and Cytoscape were used to construct the PPI network (Figure 3). Part of common DEGs was predicted to have a strong association with other genes. The size and color of the node depending on the degree, the larger the degree, the larger the node.

Especially, CCNB2, CDCA8, AURKB, BUB1B, RRM2, BIRC5, and UBE2C had the largest degree. CCNB2 is an essential component of the cell cycle regulatory machinery (Takashima et al., 2014; Li R. et al., 2019). CDCA8 is an essential regulator of mitosis and cell division (Zhang et al., 2020). AURKB participates in the regulation of alignment and segregation of chromosomes during mitosis and meiosis through association with microtubules (Ahmed et al., 2021). BUB1B encodes a kinase involved in the spindle checkpoint function (Zhang et al., 2021). RRM2 encodes one of two non-identical subunits for ribonucleotide reductase (Mazzu et al., 2020). BIRC5 encodes negative regulatory proteins that prevent apoptotic cell death (Adamopoulos et al., 2021). UBE2C is required for the destruction of mitotic cyclins and cell cycle progression (Jin et al., 2020).

To investigate the biological function of common DEGs, we used clusterProfiler to perform Functional enrichment analysis. Biological Process (BP) enrichment showed that the common DEGs were enriched in neutrophil mediated immunity, neutrophil degranulation, neutrophil activation involved in immune response, neutrophil activation and regulation of inflammatory response (Figure 4A). Cellular Component (CC) enrichment showed that the common DEGs were mainly enriched in secretory granule lumen, cytoplasmic vesicle lumen, vesicle lumen, secretory granule membrane and vacuolar membrane (Figure 4B). Molecular Function (MF) enrichment showed that the common DEGs were significantly enriched in tubulin binding, microtubule binding, carbohydrate binding, cargo receptor activity and immunoglobulin binding (Figure 4C).

KEGG pathway analysis provided a potential functional cluster of common DEGs, indicating that the common DEGs were clustered in Herpes simplex virus one infection, Human T−cell leukemia virus one infection, Cell cycle, Transcriptional misregulation in cancer and Epstein−Barr virus infection (Figure 4D).