CD36 identified by weighted gene co-expression network analysis as a hub candidate gene in lupus nephritis

Expression profile data collection

The overall procedures of our study are illustrated in the flow chart (Fig. 1). The gene expression profile of GSE104948 was selected and obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The raw data is available in the GEO database. The dataset consists of microarray-based gene expression profiles of 32 LN samples and 18 normal samples of glomerular tissues from SLE patients and living kidney transplant donors respectively. The glomerular tissues were microdissected and verified with glomerular-selective transcripts (Grayson et al., 2018). The expression values have already been log2 transformed.

Data preprocessing

The probe annotation was conducted under the Perl environment with the microarray platform file. Probes matching with multiple genes were removed, and for genes corresponding to multiple probes, the average values of probes were regarded as the expression values of the genes. A total of 11,884 genes were left after the probe annotation.

Since non-varying genes are usually regarded as background noise, we filtered the genes by variance and the top 50% (5,942 genes) with larger variance were chosen for subsequent analyses.

Weighted co-expression network construction and module division

Before WGCNA, we excluded the outlier samples by sample clustering with the hierarchical clustering method. The sample IDs and details of all samples included in our study are available in Table S1. After that, we applied WGCNA with the expression profile by using the WGCNA package (Langfelder & Horvath, 2008) in the R environment (version 3.5.3). Firstly, we calculated Pearson’s correlation for all pair-wise genes and constructed a correlation matrix. Secondly, the correlation matrix was transformed to an adjacency matrix (also known as scale-free network) with an appropriate soft-thresholding value β. A reasonable β value would emphasize strong correlations between genes and penalize weak ones. We calculated the scale-free fit index and mean connectivity of each β value from 1 to 30 respectively, and when the scale-free fit index is up to 0.85, the β value with highest mean connectivity is deemed as the most appropriate one. Then, the adjacency matrix was converted to a topological overlap matrix (TOM) so that the indirect correlations between genes are concerned. Finally, we used average linkage hierarchical clustering according to the TOM-based dissimilarity measure to classify all genes into several co-expression modules with a minimum size of 30 genes, thereby the genes with similar expression patterns were divided into the same module. After defining the first principal component of a given module as eigengene, we calculated Pearson’s correlations of the eigengenes, and merged modules whose eigengenes were highly correlated (with Pearson’s correlation higher than 0.75) into one module.

To verify the reliability of the division of modules, we plotted an adjacency heatmap of all the 5,942 genes analyzed by WGCNA. Besides, we completed a cluster analysis of module eigengenes and plotted an adjacency heatmap to find out the interactions among modules.

Identification of clinically significant modules

The clinical traits of our samples included normal and LN, we calculated the correlation between modules and traits. Modules of positive correlation with LN were considered as playing roles in the pathogenesis of the disease. On the other hand, genes in modules of positive correlation with normal trait are indispensable for maintaining normal biological functions. Thus, we extracted gene modules of highest correlation with LN and normal for subsequent studies.

Here, we introduced the definition of gene significance (GS) and module membership (MM), which represent the correlation of a given gene with clinical trait and module eigengene respectively. Genes in clinical-related modules should have high values and preferable correlations of GS and MM.

GO and KEGG pathway enrichment analyses

To explore the involved signal pathways and biological characteristics of genes in clinical-related modules, we conducted Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and visualized the top 10 significant terms respectively with the clusterProfiler R package (Yu et al., 2012). For both of GO and KEGG, enrichment terms arrived the cut-off criterion of p-value <0.01 and Benjamin-Hochberg adjusted p-value <0.01 were considered as significant ones.

Differentially expressed genes analysis

To investigate the difference of the expression profiles between normal and LN samples of genes in clinical-related modules, we applied differentially expressed genes (DEGs) analysis based on Empirical Bayes test with the limma R package (Ritchie et al., 2015). The cut-off criterion was set as follow: —log2fold change (logFC)—>1; p-value <0.01; false discovery rate (FDR) <0.001. The results were visualized with the heatmap R package (Galili et al., 2018).

Identification of hub gene

Hub gene of the LN-related module should have high connectivity with the whole module and LN trait, which may play critical roles in the molecular mechanism of LN. For identifying the hub gene related with LN, we extracted gene clusters that enriched in certain GO terms from the WGCNA network to construct sub-networks after GO enrichment analysis of the LN-related module. Then, we utilized the Cytoscape software and its plug-in cytohubba to seek out the hub gene from sub-networks (Shannon et al., 2003). After calculating the Maximal Clique Centrality (MCC) value of each gene, those with high MCC values was regarded as hub genes (Chin et al., 2014). The results were exhibited with the Cytoscape software. We then surveyed the GS value, MM value and logFC value of the selected hub gene to validate its reasonability.

Validation of hub gene with the other GEO datasets

To further verify the differential expression level of the hub gene between normal and LN tissues, we analyzed the logFC value of the hub gene with data from the other three GEO datasets (GSE32591, GSE99339 and GSE113342). The GEO IDs and details of the datasets were given in Table S1.

Validation of the clinical significance of hub gene by Nephroseq database

To assess the relationship between the expression level of the hub gene and the activity or grade of LN, we visited the Nephroseq database (http://v5.nephroseq.org), which provides unique access to datasets from the Applied Systems Biology Core at the University of Michigan, incorporating clinical data which is often difficult to collect from public sources. We then analyzed the difference of the expression level of hub gene between patients in different WHO Lupus Nephritis Class (Weening et al., 2004) based on two datasets (the dataset of Peterson Lupus Glom and the dataset of Berthier Lupus Glom) from the Nephroseq database (details available in Document S4). We performed unpaired t test for comparisons between groups and set the criterion of two-tailed value of p <0.05 as statistically significant.

Data preprocessing

After data preprocessing, 5,942 genes were selected for subsequent analyses. Sample clustering excluded the outlier samples and a total of 18 normal samples and 21 LN samples were left. The final result of sample clustering is shown in Fig. 2A, revealing satisfactory intra-group consistency and distinct difference between groups. The GEO IDs and details about the source of the samples are available in Table S1.

Weighted co-expression network construction and module division

In the current study, taking both scale-free fit index and mean connectivity as reference, the soft-thresholding was determined as 10 (Figs. 2B and 2C). Accordingly, the correlation matrix was transformed to an adjacency matrix and then converted to a topological overlap matrix. Based on average linkage hierarchical clustering and module merging, the genes were divided into six modules and were displayed with different colors (Fig. 3A), including the black, blue, brown, magenta, pink, and grey modules, containing 220, 2,881, 777, 465, 770, and 829 genes, respectively. Genes in the grey module were those couldn‘t be divided into any co-expression modules.

Figure 3B depicts the topological overlap adjacency among all the 5,942 genes analyzed by WGCNA, indicating that most genes have higher correlation with genes in the same module and lower correlation with genes in other modules, which means the division of the modules was accurate. The clustering dendrogram and adjacency heatmap of eigengene are shown in Figs. 3C and 3D, meaning that the five modules were mainly separated into two clusters.

Identification of clinically significant modules

We calculated the module-trait correlation coefficients and showed the results in Fig. 4A. The results illuminated that the blue module displayed highest correlation with LN trait (r = 0.91, p = 2e − 15), while the brown module related best with normal trait (r = 0.66, p = 4e − 06). The GS and MM value of all member genes of the blue and brown modules were shown in the scatterplots (Fig. 4B and Fig. 4C). The GS and MM value were of high correlation in the two modules (cor = 0.89, p <1e−200, and cor = 0.73, p = 3.1e−130 respectively), suggesting that the genes in the two modules were associate with their module eigengenes and clinical traits synchronously and thus suitable for further analyses and hub gene excavation. We then renamed the two modules as top LN module and top non-LN module respectively.

DEGs analysis of trait-related modules

We applied DEGs analysis for the two trait-related modules. For the top LN module, 203 DEGs were screened in LN compare with normal, including 195 up-regulated genes and 8 down-regulated ones. While, the top non-LN module contained 78 DEGs, all of which were down-regulated. The top 30 of up-regulated and down-regulated genes are displayed in Fig. 5, respectively. The logFC value, p- value, and FDR of DEGs are given in Table S2.

GO and KEGG Enrichment analyses of trait-related modules

To confirm the biological themes of genes in the trait-related modules and find the underlying biological pathways behind LN, we performed GO and KEGG enrichment analyses towards the top LN and top non-LN modules and the top 10 significant terms of GO and KEGG are exhibited in Figs. 6 and 7 respectively. The complete results of GO and KEGG enrichment analyses were given in Table S3.

For the top LN module, enriched GO-BP terms were mainly about activation of immune response, immune cells, cytokine production, and inflammation (Fig. 6A), such as “neutrophil activation” (gene count = 178, p = 7.26e−301.52E-20), “positive regulation of defense response” (gene count = 160, p = 9.14e−25), “activation of innate immune response” (gene count = 107, p = 3.45e−18). Enriched GO-MF terms were mainly about cytokine and their receptors (Fig. 6B), such as “cytokine binding” (gene count = 39, p = 2.47e−8), “cytokine receptor activity” (gene count = 34, p = 5.3e−7), “cytokine receptor binding” (gene count = 69, p = 2.98e−5). For GO-CC, enriched terms were mainly involved in membrane and endocytic vesicle (Fig. 6C), such as “vesicle lumen” (gene count = 107, p = 4.13e−15), “cytoplasmic vesicle lumen” (gene count = 106, p = 8.99e−15), “membrane region” (gene count = 95, p = 2.33e−12). The results of KEGG enrichment were similar to that of GO-BP, were mainly about immune and inflammation (Fig. 7A). The top KEGG terms included “Complement and coagulation cascades” (gene count = 31, p = 6.55e−6), “Fc gamma R-mediated phagocytosis” (gene count = 34, p = 1.90e−5), “Th1 and Th2 cell differentiation” (gene count = 33, p = 3.05e−5). Meanwhile, several well-known pathways involved in LN were also included, such as Th17 cell differentiation.

On the other hand, for the top non-LN module, the GO-BP results were mainly involved in the metabolism process of different kinds of molecule (Fig. 6D), for example, “small molecule catabolic process” (gene count = 93, p = 3.27e−41), “organic acid catabolic process” (gene count = 73, p = 3.97e−39), “carboxylic acid catabolic process” (gene count = 73, p = 3.97e−39). The top 10 terms of GO-MF and GO-CC are also displayed (Fig. 6E and Fig. 6F). Similarly, the top KEGG terms were mainly about metabolism process (Fig. 7B).

Besides, we found that the magenta module also had high correlation with LN trait. The results of enrichment analyses of magenta module are also given in Table S3.

Identification and validation of hub gene

We extracted the genes enriched in two of the top 10 significant GO-BP terms in top LN module, namely “neutrophil activation” and “positive regulation of defense response”, and constructed two sub-networks of the weighted co-expression network respectively. Co-expression pairs with top 500 weighted correlations in the sub-networks were selected for hub gene excavation. After importing the gene co-expression pairs and their weighted correlations into Cytoscape, we calculated the MCC value of genes and the most central genes of the sub-network were screened out as shown in Fig. 8. In both sub-networks, CD36 had the maximum MCC value among all, and was therefore deemed as the hub gene under the pathogenesis of LN.

The GS and MM value of CD36 were 0.772 and 0.833, revealing CD36 was significant correlated with the top LN module and LN trait. The DEGs analysis showed that the expression level of CD36 in LN was abnormally up-regulated compared with that of normal. The logFC of CD36 was 2.298, which ranked 11th among all.

Validation of hub gene with the other GEO datasets

For verifying our conclusion in a broader range, we interrogated the GEO database for more datasets about LN. We downloaded the datasets of GSE32591, GSE99339 and GSE113342, and then analyzed the differentially expressed level of CD36 between LN and normal (Fig. 9). In the four different datasets, the expression level of CD36 was consistently up-regulated in LN samples, illustrating a satisfactory reliability of the result.

Validation of hub gene by Nephroseq database

We analyzed the expression level of CD36 in glomerular tissues under different severity of LN evaluated by WHO Lupus Nephritis Class. The result indicates that CD36 was significantly up-regulated with the aggravation of LN (Fig. 10); that is to say, CD36 plays an important role in the development of LN.