Dysregulation of Signaling Pathways Due to Differentially Expressed Genes From the B-Cell Transcriptomes of Systemic Lupus Erythematosus Patients – A Bioinformatics Approach

Systemic lupus erythematosus (SLE) is an autoimmune inflammatory disorder that is clinically complex and has increased production of autoantibodies. Via emerging technologies, researchers have identified genetic variants, expression profiling of genes, animal models, and epigenetic findings that have paved the way for a better understanding of the molecular and genetic mechanisms of SLE. Our current study aimed to illustrate the essential genes and molecular pathways that are potentially involved in the pathogenesis of SLE. This study incorporates the gene expression profiling data of the microarray dataset GSE30153 from the Gene Expression Omnibus (GEO) database, and differentially expressed genes (DEGs) between the B-cell transcriptomes of SLE patients and healthy controls were screened using the GEO2R web tool. The identified DEGs were subjected to STRING analysis and Cytoscape to explore the protein–protein interaction (PPI) networks between them. The MCODE (Molecular Complex Detection) plugin of Cytoscape was used to screen the cluster subnetworks that are highly interlinked between the DEGs. Subsequently, the clustered DEGs were subjected to functional annotation with ClueGO/CluePedia to identify the significant pathways that were enriched. For integrative analysis, we used GeneGo MetacoreTM, a Cortellis Solution software, to exhibit the Gene Ontology (GO) and enriched pathways between the datasets. Our study identified 4 upregulated and 13 downregulated genes. Analysis of GO and functional enrichment using ClueGO revealed the pathways that were statistically significant, including pathways involving T-cell costimulation, lymphocyte costimulation, negative regulation of vascular permeability, and B-cell receptor signaling. The DEGs were mainly enriched in metabolic networks such as the phosphatidylinositol-3,4,5-triphosphate pathway and the carnitine pathway. Additionally, potentially enriched pathways, such as the signaling pathways induced by oxidative stress and reactive oxygen species (ROS), chemotaxis and lysophosphatidic acid signaling induced via G protein-coupled receptors (GPCRs), and the androgen receptor activation pathway, were identified from the DEGs that were mainly associated with the immune system. Four genes (EGR1, CD38, CAV1, and AKT1) were identified to be strongly associated with SLE. Our integrative analysis using a multitude of bioinformatics tools might promote an understanding of the dysregulated pathways that are associated with SLE development and progression. The four DEGs in SLE patients might shed light on the pathogenesis of SLE and might serve as potential biomarkers in early diagnosis and as therapeutic targets for SLE.

Systemic lupus erythematosus (SLE) is an autoimmune inflammatory disorder that is clinically complex and has increased production of autoantibodies. Via emerging technologies, researchers have identified genetic variants, expression profiling of genes, animal models, and epigenetic findings that have paved the way for a better understanding of the molecular and genetic mechanisms of SLE. Our current study aimed to illustrate the essential genes and molecular pathways that are potentially involved in the pathogenesis of SLE. This study incorporates the gene expression profiling data of the microarray dataset GSE30153 from the Gene Expression Omnibus (GEO) database, and differentially expressed genes (DEGs) between the B-cell transcriptomes of SLE patients and healthy controls were screened using the GEO2R web tool. The identified DEGs were subjected to STRING analysis and Cytoscape to explore the protein-protein interaction (PPI) networks between them. The MCODE (Molecular Complex Detection) plugin of Cytoscape was used to screen the cluster subnetworks that are highly interlinked between the DEGs. Subsequently, the clustered DEGs were subjected to functional annotation with ClueGO/CluePedia to identify the significant pathways that were enriched. For integrative analysis, we used GeneGo Metacore TM , a Cortellis Solution software, to exhibit the Gene Ontology (GO) and enriched pathways between the datasets. Our study identified 4 upregulated and 13 downregulated genes. Analysis of GO and functional enrichment using ClueGO revealed the pathways that were statistically significant, including pathways involving T-cell costimulation, lymphocyte costimulation, negative regulation of vascular permeability, and B-cell receptor signaling. The DEGs were mainly enriched in metabolic networks such as the phosphatidylinositol-3,4,5-triphosphate pathway and the carnitine pathway. Additionally, potentially enriched pathways, such as the signaling pathways induced by oxidative stress and reactive oxygen species (ROS), chemotaxis and lysophosphatidic acid signaling induced via G protein-coupled receptors (GPCRs), and the androgen INTRODUCTION Systemic lupus erythematosus (SLE), also known as lupus, is a rare systemic autoimmune disease that mostly affects middle-aged women, mainly of Asian, African, American, and Hispanic origin (Costa-Reis and Sullivan, 2013;Cui et al., 2013;Gurevitz et al., 2013). SLE affects an estimated 5 million people across the world, with an incidence of 1-10 per 100,000 person-years (Pons-Estel et al., 2010). SLE is characterized by a wide range of different autoantibodies, deposition of immune complexes, and immune system infiltration and inflammation within damaged organs. SLE autoantibodies invade the patient's kidneys, heart, skin, joints, and brain, leading to various typical clinical symptoms. The most common clinical symptoms of lupus are rash, arthritis, and fatigue. Severe complications of SLE lead to nephritis, anemia, neurological symptoms, and thrombocytopenia, eventually leading to severe morbidity and mortality. SLE is characterized by its clinical heterogeneity, with a wide range of clinical manifestations reflecting its complex etiopathogenesis (Tan et al., 1982). The clinical heterogeneity of SLE highlights the contribution of genetic and environmental factors to the susceptibility to the disease (Prokunina and Alarcon-Riquelme, 2004;Harley et al., 2009;Yang and Lau, 2015;Dang et al., 2016;Wang et al., 2017). To date, the reason for phenotypic variation in SLE is unknown. Understanding the molecular mechanisms behind the pathogenesis of SLE phenotypes could help in developing more efficient therapeutic approaches and preventive strategies.
With the extensive use of gene detection methods, highthroughput sequencing and extensive microarray data profiling studies on SLE have been conducted, and several differentially expressed genes (DEGs) and cellular pathways in SLE have been identified (Borrebaeck et al., 2014;Zhu et al., 2015). Nevertheless, until now, no particular gene has been recognized to act as a potential marker for the diagnosis of SLE. In addition, a large amount of data obtained from microarray technology and highthroughput sequencing have not been fully used. Ducreux et al. (2016) collected blood samples from SLE patients and healthy volunteers to identify differentially expressed genes (Ducreux et al., 2016). However, the interactions among differentially expressed genes and key genes involved in the signaling pathways of SLE remain to be elucidated. In addition, previous studies of genetic factors primarily focused on single genes; nevertheless, interactions among multiple genes may result in the multisystem invasion characteristics observed in SLE (Smith et al., 2017). Remarkably, studies have shown that disease−associated gene expression networks have a potential role in the immune response, which highlights their mechanism and therapeutic value for SLE (Deng and Tsao, 2010;Bentham et al., 2015).
Integrating and reanalyzing the data using bioinformatics methods may help in identifying gene regulatory pathways, essential genes, and their associated networks in SLE disease, which can provide new and valuable ideas for understanding the molecular mechanisms and identifying reliable diagnostic and therapeutic targets of SLE. Therefore, in this study, we first conducted a comprehensive collection of genes associated with SLE from the GEO dataset with ID GSE30153. Then, we performed a bioinformatics analysis of these genes with the MCODE (Molecular Complex Detection), GeneGo, and ClueGO tools. To further explore the pathogenesis of SLE in a more specific manner, functions and pathways identified by the modules were used to indicate the biological processes and biochemical pathways related to the immune system. Finally, the genes potentially associated with arthritis, pleurisy, and myocarditis, which are the common complications of SLE, were compared with SLE-related genes to identify common genes that participated in the development of SLE. To interpret the biological relevance of these changes in gene expression, we analyzed the microarray data via an integrated bioinformatic analysis expanding on traditional microarray analysis methods, namely, Gene Ontology (GO) and pathway analysis, thereby allowing the construction of interaction networks that might identify novel prognostic markers and therapeutic targets.

Acquisition of Array Data and Processing
Gene expression profiling data from microarray array analysis of the GSE30153 dataset were downloaded from the NCBI GEO database (Gene Expression Omnibus database) 1 . The database accommodates gene expression datasets from a variety of experiments, such as DNA-seq, ChIPs, RNA-seq, microarray, and high-throughput hybridization array (Edgar et al., 2002;Barrett et al., 2013). GSE30153 contains 26 samples, including 17 patients with SLE and 9 healthy controls of human sorted B-cells obtained by using the platform GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array (Garaud et al., 2011). The downloaded gene expression profiling data are freely available in the public database, and there were no human or animal experiments conducted by any of the authors in this study.

Preprocessing of Data and DEG Identification
Using the robust multiarray standard model, the initial information from the dataset was subjected to quantile normalization, background correction, and log transition (Irizarry et al., 2003). Preprocessing included changing to gene symbols from probe IDs using the Gene ID converter from Entrez (Alibés et al., 2007). The statistical online tool GEO2R uses the R/Bioconductor, and limma package v3.26.8 was used to screen the raw gene expression data (Smyth, 2005;Barrett et al., 2013;Ritchie et al., 2015). We performed a Benjamini-Hochberg test (to determine the false discovery rate) and T-tests to compute the false discovery rate (FDR) and p-values to identify the DEGs between SLE patients and healthy control human sorted B-cells (Benjamini and Hochberg, 1995;Aubert et al., 2004). We set the primary criteria of | log (2 fold change) | > 1 and p < 0.05 to obtain significant DEGs from the dataset, whereas cutoffs of log2FC ≥ 1 and log2FC ≤ −1 were used to denote upregulated and downregulated DEGs, respectively. For high-throughput sequencing, a logarithm to base 2 is widely used and in the initial scaling, the doubling is equivalent to a log2FC of 1 (Love et al., 2014). A volcano plot was constructed using a web-based tool 2 . The resulting DEGs were used for further analysis.

Constructing PPI Networks
To assess the relationships between the DEGs from the GSE30153 dataset, we constructed a protein-protein interaction (PPI) network by using Search Tool for the Retrieval of Interacting Genes (STRING v11.0) 3 (Szklarczyk et al., 2017(Szklarczyk et al., , 2019. The cutoff criterion was set to a high confident interaction score of ≥0.7 to eliminate inconsistent PPIs from the dataset. We then incorporated the results from the STRING database into Cytoscape software (v3.7.2) 4 to envisage the PPIs within the statistically relevant DEGs (Shannon et al., 2003). The MCODE plugin from Cytoscape was utilized to identify the interconnected regions or clusters from the PPI network. The cluster finding parameters were adopted, such as a degree cutoff of 2, a node score cutoff of 0.2, a kappa score (K-core) of 5, and a max depth of 100, which limits the cluster size for coexpressing networks (Bader and Hogue, 2003). The top clusters from MCODE were subjected to ClueGO v2.5.5/CluePedia v1.5.5 analysis to obtain comprehensive GO and pathway results from the PPI network. ClueGO combines GO and pathway analyses from KEGG and BioCarta and provides a fundamentally structured GO or pathway network from the PPI network (Bindea et al., 2009).

Metacore GeneGo Analysis of DEGs
Metacore, a Cortellis Solution software (Clarivate Analytics, London, United Kingdom) 5 , was used to perform curated pathway enrichment analysis and GO analysis. GeneGo facilitates the rapid assessment of metabolic pathways, protein biological networks, and pathway maps from high-throughput experimental data (MetaCoreLogin | Clarivate Analytics). Based on a significance threshold of p < 0.05, a pictorial representation of the molecular interactions of DEGs from the study groups is generated. Determination of a hypergeometric p-value enables the estimation of the chance that an intersection between DEGs and ontological elements is random. An FDR < 0.05 was used as a criterion to calculate if statistically significant DEGs constituted a processor pathway.

Identification of DEGs From the Dataset
Our study contained the gene expression profiles of the GSE30153 dataset from the GEO database, which were submitted by Garaud et al. (2011) based on analysis with the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) (Garaud et al., 2011). The dataset encompasses 26 samples, including 17 patients with SLE and 9 healthy controls ( Table 1). By utilizing the GEO2R online tool, we obtained the differentially expressed genes (DEGs) from the GSE30153 dataset by comparing the SLE samples with control samples. By calculating p-values and | log2FC | values, the top 250 DEGs were identified. A volcano plot was constructed using the Rstudio web server ShinyVolcanoPlot to identify DEGs by comparing the SLE and control groups from the dataset. The volcano plot in Figure 1 depicts all the DEGs with a log2FC against the -log10 (p-value) between the two groups. With cutoffs of p < 0.05 and log2FC ≥ 1.0 or ≤−1, we found 4 and 13 genes that were upregulated and downregulated, respectively, between the two groups ( Table 2). The genes that were differentially expressed between the two groups are shown in Supplementary Table S1.

Screening of Module and Construction of Interlinking PPI Network
To assess the protein-protein connections among the DEGs, we used the STRING tool to compute the protein interactions and plotted them using Cytoscape v3.7.2. Figure 2 depicts the PPI network with 103 nodes and 201 edges. The DEGs are represented as nodes, and the edges are interactions between the DEGs. A combined node score of >0.4 was considered to be significant. MCODE plugin v1.5.1 from Cytoscape was utilized to identify the densely interlinked regions within the protein network. As a result, we obtained the top two significant clusters from the DEGs protein network with MCODE scores of 5.043 and 3.625. A graphical representation of these clusters is shown in Figures 3A,B. The subnetworks, scores, number of nodes and edges, and node IDs are tabulated in Table 3.  Frontiers in Bioengineering and Biotechnology | www.frontiersin.org FIGURE 2 | The network demonstrates the protein-protein interaction between the DEGs identified from GSE30153 using Cytoscape. The nodes represented as ellipse (robin's blue) and edges as lines (gray).

Enrichment Analysis by ClueGO
The top two subnetworks from MCODE were used as an input for analyzing the functional enrichment of PPI subnetworks using the ClueGO/CluePedia plugin from Cytoscape. In a biologically clustered subnetwork, ClueGO helps to visualize the biological terms of broad gene clusters. The subnetwork enrichment analyses of MCODE cluster 1 and cluster 2 are depicted in  complications ( Figure 4A). The DEGs from cluster 2 were mainly enriched in the regulation of the transcription involved in the G1/S transition of the mitotic cell cycle (GO: 0000083), the negative regulation of the G0 and G1 transitions (GO: 0070317), and the p53 signaling pathway (KEGG: 04115) ( Figure 4B).
The pathways that were activated in the enrichment analysis were highly related to B-cell pathophysiology, resulting in events associated with the immune system, vasculopathy, and kidney.

Metacore TM GeneGo for Enrichment Analysis of DEGs
Further functional enrichment analysis was carried out using Metacore TM GeneGo software from Clarivate Analytics to comprehensively dissect the pathways associated with the DEGs. Using the functional ontology feature in GeneGo, the IDs of potential genes that were involved in the target pathways were identified. Based on hypergeometric p-values, the probability that the intersection of a gene set and associated ontological objects was random was evaluated. A decreased p-value indicated that the entity would be more significant to the DEGs, suggesting a better score. The functional enrichment analysis of the DEGs defined the top 10 metabolic networks, and canonical pathway maps are depicted in Figures 5A,B. For each classification, the significant statistical data rely on a low p-value. The pathway maps with the lowest p-value are shown in Figures 6A-C. These are the top-scoring signaling pathways based on the gene enrichment distribution, which emphasizes that the DEGs from human sorted B-cells are triggered via oxidative stress and ROS-induced cellular signaling (Figure 6A), chemotaxis and lysophosphatidic acid signaling via GPCRs (Figure 6B), and androgen receptor activation and downstream signaling in prostate cancer ( Figure 6C). The well-distinguished proteins and complexes of proteins are shown as specific symbols 6 ; all experimental data are displayed and have corresponding thermometer-like symbols on all the maps. The upregulated genes are indicated by a red thermometer facing upwards.  cluster 2 (B). The plugin provides a combined enrichment analysis of clusters, including the GO biological process, molecular function, and pathway from KEGG. The GO term/pathway network connectivity defined by edges and functional clusters on genes shared between terms (kappa score ≥ 0.4) and displaying pathways only with p≤ 0.05. The size of the node indicates the p-value. The color code of nodes represents the functional group that they belong to. The most important functional terms specify the pathway names within each class are indicated in bold colored characters. (A) The network enrichment analysis of cluster 1. Each node constitutes a precise term for cluster 1; (B) The network enrichment analysis of cluster 2. Each node constitutes a precise term for cluster 2.
FIGURE 5 | The top 10 metabolic networks and pathway maps were annotated using GeneGo enrichment analysis for the genes that are differentially expressed from SLE patients vs. healthy controls, respectively. (A) The content of these metabolic networks was annotated and defined by GeneGo Cortellis Solution software. Each process represents a pre-set network of protein interactions characteristic for the process, and sorting was performed for the metabolic networks that are statistically significant. (B) The pathway maps (canonical) of GeneGo display a series of signals and metabolic charts that cover human in a structured manner. The significant expression of a gene/protein represented in histogram height.

DISCUSSION
DNA microarrays and next generation sequencing (NGS) approaches are high-throughput technologies that have resulted in the emergence of new biomedical discoveries. Data from microarray and gene expression profiles have enabled a deeper understanding of the intrinsic molecular pathways of complex mechanisms of biological systems and their responses (Russo et al., 2003;Babu, 2004;Perez-Diez et al., 2013;Kumar et al., 2019). It is therefore highly relevant to examine the peripheral B-cell transcriptomes of SLE patients and healthy controls to determine genes that are differentially regulated and their target pathways. Our current study extracted DEGs from 17 SLE patients and 9 healthy controls from the GEO database (GSE30153) (Garaud et al., 2011). The top 250 DEGs were identified, including 4 upregulated and 13 downregulated genes from the groups through bioinformatics strategies ( Table 2 and  Supplementary Table S1). These identified DEGs were subjected to ClueGO and GeneGo Metacore TM analysis for GO and pathway annotation, and constructed the interacting networks of PPI and used for cluster analysis. In the network, the nodes were considered proteins, and the edges were their interactions. Using network topology features, the PPI network can be analyzed to distinguish the core proteins that are involved in the pathways (Barabási and Oltvai, 2004;Ideker and Sharan, 2008;Keskin et al., 2016;Kumar et al., 2019). The identified DEGs from the present study were analyzed with STRING to exploit the complex interactions between the DEGs via text mining, evidence from experiments, and repositories (Figure 2). We performed module screening of the PPI networks using the MCODE plugin from Cytoscape. As a result, we obtained significant clusters that are densely interlinked regions in the PPI network (Figures 3A,B). Screening of these clusters from the network might help to identify the essential genes that are involved in the pathogenesis and progression of SLE. The obtained clusters mostly contained protein complexes or proteins present in the pathways in the PPI network, and cluster visualization is essential for comprehending the properties of the network functionally and systematically (Krogan et al., 2006;Rahman et al., 2013). Furthermore, to identify the functional enrichment of these subnetworks from MCODE, we implemented the ClueGO plugin for analysis. This revealed that the DEGs were enriched in most essential pathways, which are highly associated with the immune system. The GO and KEGG enrichment analyses of the DEGs from cluster 1 showed that they were mostly enriched in T-cell costimulation, lymphocyte costimulation, negative regulation of vascular permeability, the metaphase/anaphase transition of the mitotic cell cycle, regulation of the transcription involved in the G1/S transition of mitotic cell cycle, the hematopoietic cell lineage, the B-cell receptor signaling pathway, the ErbB signaling pathway, the AGE-RAGE signaling pathway in diabetic complications, and pancreatic cancer. Interestingly, the costimulation of T-cell and lymphocyte receptors is recognized to be important in SLE pathogenesis by enabling communicating with B-cells for the production of autoantibodies (Shlomchik et al., 2001;Mak and Kow, 2014). In SLE, negative regulation of vascular permeability may be induced by different mechanisms; the dysregulated genes from the cluster 1 subnetwork might lead to endothelial cell damage and vasculopathy (Favero et al., 2014;Lee et al., 2019). The differential cell signaling results in the recruitment of various proteins and inappropriate activation of B-cells (Zhou et al., 2009;Comte et al., 2015). Oxidative stress is common in inflammatory disorders and results in the increased production of reactive carbonyl groups that are partially converted to AGEs, and the DEGs in the AGE-RAGE signaling pathway might also be involved in the accumulation of AGEs in SLE patients and lead to diabetic complications (de Leeuw et al., 2007;Li et al., 2007;Kurien and Scofield, 2008;Nienhuis et al., 2008). Interestingly, our enrichment analysis found that the identified differential expression of the genes (AKT1, VEGFA, CDK6, and MAPK9) that were involved in the risk of developing pancreatic cancer in SLE patients was due to chronic inflammation, suggesting that these genes might be involved in the pathogenesis of SLE. Our findings are therefore consistent with the roles of genes that are differentially expressed in SLE-causing pathways ( Figure 4A). The enrichment analysis of the cluster 2 subnetwork showed that the DEGs were mostly enriched in the regulation of transcription involved in the G1/S transition of the mitotic cell cycle, the negative regulation of the G0 and G1 transitions, and the p53 signaling pathway. It has been reported that the proliferation of T-cells is followed by lowered levels of cyclin-dependent kinase (CDK) inhibitors, and alterations in the expression of CDKs in the G0/G1 phase were seen in the lymphocytes of SLE patients (Yamauchi and Bloom, 1997;Tang et al., 2009). The DEGs involved in the cluster 2 subnetwork might negatively regulate these pathways. Alterations in cyclin-CDK complex behavior and cyclin-dependent kinase inhibitors (CDKIs) have been reported to alter the proliferation of T-cells, oxidative stress, and immune responses (Santiago-Raber et al., 2001;Tang et al., 2009). p53 signaling is essential for various cellular mechanisms, and defects in this signaling pathway are associated with SLE development. Considerably elevated levels of p53 protein are found in SLE patients with active inflammatory disorders (Miret et al., 2003;Veeranki and Choubey, 2010). Apoptosis dysregulation appears to be another cause of SLE pathogenesis because the possible sources of autoantigens are cell debris from apoptosis in SLE, and excessive cellular senescence of the immune cells, especially T-cells, was reported in SLE patients with peripheral blood mononuclear cells (PBMCs) and skin lesions (Colonna et al., 2014;Sáenz-Corral et al., 2015). Thus, our identified DEGs (RRM2, APC, CHEK1, E2F6, TYMS, E2F7, and CDK6) from the cluster 2 subnetwork are highly related to and consistent with the members of the signaling pathways associated with the immune system, apoptosis, the cell cycle, and vasculopathy.
To clearly define the interactions between the proteins and signaling pathways examined from the interpretation of STRING, Cytoscape, MCODE, and ClueGO analyses, we implemented the GeneGo Metacore software, which incorporates extensive data on metabolic signaling pathways and their regulatory mechanisms and contains accurately complied networks of biological systems. By utilizing the GeneGo Metacore software, we obtained a detailed description of the DEGs that participate in SLE pathogenesis based on the determined p-values. Among the top 10 metabolic networks, the phosphatidylinositol-3,4,5-triphosphate pathway, O-hexadecanoyl-(L)-carnitine pathway, 1,2-didocosapentaenoylsn-glycerol 3-phosphate pathway, and 1-linoleoyl-glycerol-3phosphate pathway were profoundly enriched and significant in the SLE DEGs ( Figure 5A). The increased activity of phosphatidylinositol-3,4,5-triphosphate stimulates essential cell signaling pathways such as the pathways involved in cell division, survival, and the rapid increase in T-lymphocytes in SLE (Comte et al., 2015). PI3K (phosphatidylinositol 3-kinase) is a protein kinase that phosphorylates phosphatidylinositol 4,5-phosphate to regulate the signaling of T-lymphocytes; an increased amount of PI3K was also observed in an animal model of lupus (Liu et al., 1998;Grolleau et al., 2000;Niculescu et al., 2003;Joseph et al., 2014). Modification of the carnitine signaling pathway results in various organ failures by producing effective responses to pathogens (Famularo and De Simone, 1995;Famularo et al., 2004). Thus, the DEGs involved in the O-hexadecanoyl-(L)-carnitine pathway might lead to increased immune responses. In addition, the top three pathways associated with the DEGs of sorted B-cells from SLE patients were mostly enriched in oxidative stressand ROS-induced cellular signaling (Figure 6A), chemotaxis and lysophosphatidic acid signaling via GPCRs (Figure 6B), and androgen receptor activation and downstream signaling in prostate cancer ( Figure 6C). Recent findings have shown that oxidative stress and ROS induce molecular alterations that have adverse effects in SLE patients (Choi et al., 2016;Tsokos et al., 2016;Lightfoot et al., 2017). Elevated oxidative stress in SLE patients leads to the accumulation of higher amounts of oxidative lipoproteins, which are harmful in zebrafish models and cause additional oxidative damage to the system (Chung et al., 2007;Park et al., 2016;Lightfoot et al., 2017). Interestingly, our study identified the EGR1 gene as downregulated in the SLE patients in comparison to controls, and it also plays a role in ROS signaling. This clearly indicates that EGR1 might be required to maintain the oxidative stress and ROS signaling pathways.
Moreover, the DEGs involved in the oxidative stress signaling pathway might contribute to peripheral neuropathy, damage to blood vessels, and cardiovascular events, which are the prominent clinical conditions found in SLE patients. Chemotaxis and lysophosphatidic acid (LPA) signaling are essential pathways in autoimmune inflammatory disorders, and GPCRs are responsible for regulating immune cells via LPA receptors (Yang et al., 2005;Skoura and Hla, 2009). G2A gene knockout resulted in the hyperresponsiveness of T-cells to T-cell receptor stimulation, manifesting as an increased proliferation of T-cells, which may promote inflammatory phenotypes in G2A-deficient mice (Le et al., 2001). Various studies have suggested that LPA plays a vital role in atherosclerosis progression and development by promoting neutrophil and monocyte adherence and enhancing inflammatory events (Siess et al., 1999;Smyth et al., 2008;Skoura and Hla, 2009). The androgen receptor (AR) is a transcription factor that is activated by a ligand and is essential for cells targeted by the androgen response (Robeva et al., 2013;Gubbels Bupp and Jorgensen, 2018). AR also regulates immune function in SLE via transcriptional regulation of various genes. Our study identified the transcription factor AR, which positively regulates the c-Myc, SCAP, prosaposin, and KLF5 genes, which are responsible for inflammatory responses, and promotes tumor growth factors and cytokine signaling when activated ( Figure 6C). The enriched terms from ClueGO modules and the GeneGo-identified terms correlated well in this study and validate the significance of the findings from the pathway maps. The combined results from these two enrichment analyses suggest that B-cells from SLE patients and B-cells from healthy controls undergo differential gene expression associated with positive regulation of kidney development, the hematopoietic cell lineage, positive regulation of vasoconstriction, T-cell costimulation, and regulation of the transcription involved in the G1/S transition of the mitotic cell cycle.
In addition to the interaction analysis, we carried out interrelation analysis for the essential genes to determine the relationships between the genes, which implicitly or explicitly interacted with each other. Interestingly, the identified genes indirectly communicated with each other via molecular signaling pathways related to mTOR signaling, apoptosis, PI3K-Akt signaling, the hematopoietic cell lineage, positive regulation of vasoconstriction, signaling by receptor tyrosine kinases, AGE-RAGE signaling, and lymphocyte and T-cell costimulation (Figure 7). EGR1 and AKT1 are directly involved in oxidative stress via ROS and AGE-RAGE signaling, whereas CAV1 is directly involved in tyrosine kinase receptor signaling and lymphocyte and T-cell costimulation. CD38 is directly associated with the hematopoietic cell lineage and positive regulation of vasoconstriction. Overall, the dysregulation of the indicated pathways in SLE patients is a result of differential gene expression. The essential genes are differentially expressed between cells from patients with SLE and cells from healthy controls and are present in important signaling cascades, which could be a crucial factor for SLE development.

CONCLUSION
Taken together, the results of our comprehensive bioinformatics analysis showed that the DEGs identified between sorted B-cells from patients with SLE sorted B-cells from controls could play a significant role in the growth, progression, and development of SLE. This study identified 4 upregulated and 13 downregulated genes, including essential genes (EGR1, CD38, CAV1, and AKT1), from the pathway enrichment analysis. Indeed, the identified pathways from the enrichment analysis were strongly related to the immune system, vasculopathy, cardiovascular functions, and inflammatory responses, which are processes that can lead to the development of SLE. The broad understanding of SLE pathophysiology from this study will allow us to identify and develop therapies targeting SLE and contribute to personalized treatment strategies. Collectively, the study findings could aid in enhancing our understanding of the fundamental molecular processes of SLE and provide possible strategies for early diagnosis in SLE; in addition, combinatorial therapeutic strategies using oxidative stress and ROS cellular signaling and lysophosphatidic acid signaling via GPCRs might have symbiotic effects on the molecular events in SLE.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE30153.

AUTHOR CONTRIBUTIONS
SU, DT, HZ, and CG were involved in the design of the study and the acquisition, analysis, and interpretation of the data. SU, DT, CG, SY, NY, MS, and HZ were involved in the interpretation of the data and drafting the manuscript. CG, RS, and HZ supervised the entire study and were involved in study design, the acquisition, analysis, and interpretation of the data, and drafting the manuscript. The manuscript was reviewed and approved by all the authors.