Assessment of the genes and molecular mechanisms of B cells activation through systems biology approaches

2.1 Microarray dataset analysis

The expression profile of GSE84948 microarray dataset deposited by Rosenberg et al. was extracted from the Gene Expression Omnibus (GEO) NCBI database. The GSE84948 dataset explores gene expression profiles of B cells cultured with Th1 and Th2 cells. The quality control assessment of microarray data was calculated by unsupervised principle component analysis (PCA) method using ggplot2 package of R software [4]. GEO2R web tool of GEO was used to identify differentially expressed (DE) genes. Controls (naive B cells) were compared with Be1 and Be2 cells, respectively. The comparison of the groups carried out using student’s t-test with Benjamini–Hochberg false discovery rate p-value correction. Genes were considered as differentially expressed, had adjusted p-value less than 0.05.

2.2 Protein-protein interaction networks construction

Using CluePedia plugin version 2.1.7 [5] of Cytoscape software version 3.7.1 [6], the protein-protein interaction (PPI) networks were constructed with discriminated DE genes. For retrieving interactions, STRING database with confidence cutoff 0.80 was applied [7]. Interactions between genes were based on activation, binding, post-translational modification, and inhibition. The critical nodes with the highest score of graph theory parameters such as degree, betweenness centrality, and closeness centrality were determined by the NetworkAnalyzer tooll of Cytoscape [8]. Using Molecular Complex Detection (MCODE) plugin of Cytoscape, clustering analysis of networks were performed according to the default settings [9].

ClueGO plugin version 2.1.7 [10] of Cytoscape was used for Gene ontology (GO) analysis with genes (nodes) each network. GO “Biological Process” and adjusted p-value ⩽ 0.05 were chosen for retrieving pathways [11]. For p-value correction Bonferroni step down was considered.

3.1 Differentially expressed genes were determined after quality assessment

In this study, we re-analyzed the GSE84948 dataset, the expression profiles contain naive B cells and B cells activated with Th1 and Th2 (Be1 and Be2 cells). In the quality analysis step, the PCA plot showed that groups were segregated based on their state, indicating the satisfactory quality of this dataset (Fig. 1). The comparison of naive B cells with Be1 and Be2 cells identified 8742 and 8748 DE genes with adjusted P-value ⩽ 0.05, respectively (supplementary Table 1). 7767 overlapping genes between groups were ignored for future analysis.

OPEN IN VIEWER

The PPI networks constructed with 975 and 981 genes from Be1 and Be2 cells analyses, respectively, achieved complex networks with 469 and 460 significant nodes. The topology of the networks were analyzed to identify the crucial genes. Nodes with high value of the degree, betweenness centrality, and closeness centrality parameters were illustrated as critical genes. Jak3, Actrt3, and Pik3cb genes were determined as central genes in Be1 network. Prkx, Smarca4, and Jak2 genes are the key genes in Be2 PPI network. The 15% of top genes were shown in Table 1. These genes were also detected in modules obtained from cluster analysis.

3.3 Gene ontology (GO) and enrichment analysis were performed

Using the PPI networks genes, Gene ontology (GO) analysis was carried out and we reach to processes that strongly related and connected to B cell development. Response to cytokine, lymphocyte activation and differentiation are the same functions in both analyses. MAPK cascade, ERK1 and ERK2 cascade, and NIK/NF-kappaB signaling are pathways obtained from Be1 GO analysis. STAT cascade, response to interleukin-1, and response to molecule of bacterial origin were found from Be2 assay. Table 2 shows the results of GO analysis.