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

Abstract

Two classes of T helper lymphocytes, Th1 and Th2, have different roles in B cell activation based on specific cytokines. To understand the difference of molecular mechanisms of B cell activation, the microarray dataset of B cells co-cultured with type 1 and 2 T helper, Be1 and Be2, were investigated. After quality assessment, using the GEO2R tool, the GSE84948 dataset was re-analyzed. Genes with adjusted $p$ -value $\leqslant$ 0.05 were assumed as differentially expressed (DE). The protein-protein interaction (PPI) networks were constructed using CluePedia plugin of Cytoscape, and analyzed by NetworkAnalyzer tool and MCODE plugin. Using ClueGO plugin of Cytoscape software, gene ontology (GO) analysis was performed. The comparison of Be1 and Be2 cells with naive B cells revealed 8742 and 8748 DE genes, respectively. The topology analysis of PPI networks predicted central genes. Among these, Jak3, Actrt3, and Pik3cb genes were determined as crucial genes in Be1 network. Prkx, Smarca4, and Jak2 genes were defined as critical genes in Be2 PPI network. GO analysis with PPI networks genes resulted in the promotion of immune system activation. In conclusion, we explored holistic methods for molecular assay of the difference between B cell activation mechanisms with Th1 and Th2.

Keywords

T-Lymphocytes helper-Inducer cytokines microarray protein interaction network gene ontology

Table 1
Central genes were determined. The top 15% genes with the highest scores of degrees, betweenness centrality (BC), and closeness centrality (CC) in PPI networks were obtained

Be1						Be2
GOTerm	Degree	GOTerm	BC	GOTerm	CC	GOTerm	Degree	GOTerm	BC	GOTerm	CC
Jak3	151	Actrt3	0.09	Pik3cb	0.43	Prkx	156	Smarca4	0.12	Jak2	0.45
Matk	130	Akt1	0.08	Mtor	0.42	Jak2	143	Jun	0.09	Itk	0.44
Rhobtb1	112	Jak3	0.07	Jak3	0.41	Itk	134	Pik3cd	0.09	Pik3cd	0.43
Rap1b	112	Mtor	0.07	Rhobtb1	0.41	Pik3cd	106	Prkx	0.08	Smarca4	0.43
Akt1	111	Ppp3cc	0.06	Akt1	0.41	Cdk2	102	Jak2	0.07	Actr1b	0.42
Rhobtb3	111	Nfyb	0.06	Rhobtb3	0.41	Rps6ka1	97	Dvl2	0.06	Rap1a	0.42
Rhebl1	109	Rhobtb1	0.04	Prkcq	0.40	Rnd2	91	Itk	0.06	Cdk2	0.42
Pik3cb	108	Pkp2	0.04	Actrt3	0.40	Smarca4	90	Actr1b	0.05	Rnd2	0.42
Prkch	94	Rhobtb3	0.04	Prkch	0.40	Phlpp1	87	Cdk2	0.05	Rap2a	0.42
Prkcq	92	Matk	0.03	Matk	0.40	Rap1a	85	Phlpp1	0.05	Rasd1	0.41
Mtor	88	Ddx39b	0.03	Rap1b	0.40	Rap2a	82	Rps6ka1	0.04	Phlpp1	0.41
Ppp3cc	87	Prkch	0.03	Rhebl1	0.40	Rasd1	81	Itgb3	0.04	Jun	0.41
Nfyb	86	Rap1b	0.03	Ppp3cc	0.39	Ptpn6	75	Ptpn6	0.04	Rps6ka1	0.41
Rps6kb2	85	Cdk5	0.03	Rps6kb2	0.39	Itgb3	68	Gnai3	0.04	Prkx	0.41
Rps6ka4	79	Ifng	0.03	Cdk5	0.39	Actr1b	66	Polr2e	0.03	Itgb3	0.40
Pim2	74	Rhebl1	0.03	Chd1l	0.38	Jun	63	Ppp1cb	0.03	Gnai3	0.39
Araf	72	Rxra	0.03	Rps6ka4	0.38	Gnai3	61	Rap1a	0.03	Ptpn6	0.39
Actrt3	69	Chd1l	0.03	Araf	0.38	Sfn	60	Edn1	0.03	Gnaz	0.39
Grap	67	Gnb4	0.03	Grap	0.38	Pim3	58	Ubqln4	0.03	Lpxn	0.39
Arf1	62	Uty	0.03	Chd1	0.38	Tyro3	50	Rnd2	0.03	Edn1	0.39
Arl5b	60	Prkcq	0.02	Dnm3	0.38	Ern1	48	Rap2a	0.02	Nfkb1	0.38
Cdk5	55	Setd1a	0.02	Gnb4	0.37	Gnaz	47	Atp6v1f	0.02	Sfn	0.38

1. Introduction

Two major subtypes of T helper lymphocytes are known as type 1 and 2, (Th1) and (Th2), respectively [1]. Th1 cells respond against intracellular parasites, whereas Th2 cells cause immune responses against extracellular organisms. They secrete specific cytokines in response to antigenic stimulation. Th1 cells produce interleukin (IL)-2 and interferon (IFN)-g, while Th2 cells produce IL-4, IL-5, IL-6, IL-10, and IL-13 [2, 3]. Our objective was to clearify the difference of B cell activation mechanisms into producing antibody secreting cells in the presence of Th1 and Th2. For this purpose, we used system biology approaches to reach a comprehensive understanding of each process. We re-analyzed the microarray expression profiling of activated B cells co-cultured with Th1 and Th2. Be1 cells (B cells cultured with Th1 cells) and Be2 cells (B cells cultured with Th2 cells) were compared with naive B cells. Beside functional analysis, with protein-protein interaction (PPI) network analysis, we identified novel genes related to B cell differentiation.

2. Methods

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].

Table 2
Gene ontology (GO) analysis was carried out with genes (nodes) in PPI networks. Results shown the terms which are known to be related with B cell activation. Number of genes (Nr. Genes) for each term mentioned

Be1		Be2
GOTerm	Nr. genes	GOTerm	Nr. genes
Anatomical structure morphogenesis	116	Establishment of localization in cell	79
Cellular response to organic substance	106	Cellular macromolecule localization	67
Cellular response to stress	83	Cytoskeleton organization	62
Tissue development	83	Peptidyl-amino acid modification	62
Establishment of protein localization	78	Response to drug	53
Peptide transport	75	Response to cytokine	52
Cell motility	70	Immune system development	51
Cell migration	66	Leukocyte activation	50
Cellular response to endogenous stimulus	62	Hematopoietic or lymphoid organ development	48
Peptidyl-amino acid modification	55	Response to lipid	47
Cell morphogenesis	53	Lymphocyte activation	42
Leukocyte activation	53	Cell-cell adhesion	41
MAPK cascade	52	Inflammatory response	37
Cellular response to oxygen-containing compound	50	Supramolecular fiber organization	37
Secretion by cell	46	Apoptotic signaling pathway	33
Lymphocyte activation	45	Leukocyte differentiation	33
Response to cytokine	44	Protein localization to membrane	30
Innate immune response	42	Chemotaxis	29
Apoptotic signaling pathway	40	Adaptive immune response	28
Cellular response to nitrogen compound	35	Homeostasis of number of cells	27
Inflammatory response	35	Leukocyte cell-cell adhesion	26
Leukocyte differentiation	34	Lymphocyte proliferation	26
Response to growth factor	33	Lymphocyte differentiation	25
Adaptive immune response	31	Response to molecule of bacterial origin	24
Lymphocyte differentiation	24	Cellular component disassembly	22
ERK1 and ERK2 cascade	23	Leukocyte migration	21
Cell-substrate adhesion	21	Protein polymerization	19
Endosomal transport	14	Cellular response to biotic stimulus	16
DNA-templated transcription, initiation	13	Myeloid cell homeostasis	16
NIK/NF-kappaB signaling	13	Leukocyte apoptotic process	15
Androgen receptor signaling pathway	8	STAT cascade	15
Collateral sprouting	7	Interferon-gamma production	12
		Response to interleukin-1	12
		CD4-positive, alpha-beta T cell activation	11
		Leukocyte homeostasis	11
		Response to protozoan	7
		T-helper cell lineage commitment	5

2.3 Gene ontology (GO) analysis

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 $\leqslant$ 0.05 were chosen for retrieving pathways [11]. For p-value correction Bonferroni step down was considered.

3. Results

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 $\leqslant$ 0.05, respectively (supplementary Table 1). 7767 overlapping genes between groups were ignored for future analysis.

Figure 1.

The quality of the microarray dataset is satisfying. The Principle component analysis result of GSE84948 dataset which consist of the expression naïve B cells, Be1 cells, Be2 cells is acceptable.

3.2 Central genes in protein-protein interaction networks were identified

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.

4. Discussion

T helpers play important roles in immune responses to intracellular and extracellular pathogens [3]. The roles of Th1 and Th2 during the B cell activation are not exactly defined. Here, with the new high-throughput methodology, we discovered that main genes have been associated with B cell activation with Th1 and Th2 cells, Be1 and Be2, respectively. Furthermore, with gene ontology (GO), the related biological processes were shown. The map interaction network assay and function analysis of Be1 and Be2 cells do not show completely different and exhibit some distinctions.

The PPI network analysis investigates the genes interactions. To determine the critical nodes in the complex networks, we employed a combination of different measures of centrality. Some nodes such as Jak3 and Prkx have high degree, in Be1 and Be2 networks, respectively, so they have many connections and are vital for the surveillance of the networks. Betweenness centrality measures the number of shortest paths going through a node. Nodes with high between ness centrality such as Actrt3 and Smarca4 in Be1 and Be2 networks, respectively, are shortcuts of the networks. In addition, nodes with highest closeness centrality such as Pik3cb and Jak2 in Be1 and Be2 networks, respectively, are physically nearest genes to all nodes [8, 12]. Also, modules were explored to identify functional genes. Modules are high density regions in the network. All of Actrt3, Akt1, Cdk5, Jak3, Matk, Mtor, Ppp3cc, Prkch, Prkcq, Rap1b, Rhebl1, Rhobtb1, Rhobtb3, and Rps6kb2 genes have the highest scores of centrality in the Be1 PPI network. Actr1b, Cdk2, Edn1, Ern1, Gnai3, Itgb3, Itk, Jak2, Jun, Phlpp1, Pik3cd, Prkx, Ptpn6, Rap1a, Rap2a, Rnd2, Rps6ka1, and Smarca4 genes were determined as critical genes in Be2 PPI network. Moreover, these genes were found in modules. Rps6kb2 and Edn1 have the highest rank of all centrality parameters, also detected as seed genes in cluster analysis, in Be1 and Be2 PPI networks, respectively. For Be1 and Be2, 32 and 37 significant functions were obtained from 469 and 460 nodes, respectively, most of them involved in promoting immune system.

5. Conclusions

In this study, using computational tools we generate a systematic view of the B cell activation in the presence of Th1 and Th2 to determine the genes basic molecular mechanism of B cell-activating functions. The role of these genes and functions remains to be confirmed in future experimental studies.

Footnotes

Acknowledgments

This study was supported by the Isfahan University of Medical Science.

Conflict of interest

The author declares no competing financial interests.

Supplementary data

The supplementary files are available to download from https://dx-doi-org.web.bisu.edu.cn/10.3233/HAB-190393.

References

Mosmann

T.R.

Cherwinski

Bond

M.W.

Giedlin

M.A.

and Coffman

R.L.

, Two types of murine helper T cell clone. I. Definition according to profiles of lymphokine activities and secreted proteins, J Immunol 136 (1986), 2348–2357.

Mosmann

T.R.

and Coffman

R.L.

, TH1 and TH2 cells: different patterns of lymphokine secretion lead to different functional properties, Annu Rev Immunol 7 (1989), 145–173.

O’Garra

and Arai

, The molecular basis of T helper 1 and T helper 2 cell differentiation, Trends Cell Biol 10 (2000), 542–550.

Yeung

K.Y.

and Ruzzo

W.L.

, Principal component analysis for clustering gene expression data, Bioinformatics 17 (2001), 763–774.

Bindea

Galon

and Mlecnik

, CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data, Bioinformatics 29 (2013), 661–663.

Shannon

Markiel

Ozier

Baliga

N.S.

Wang

J.T.

Ramage

Amin

Schwikowski

and Ideker

, Cytoscape: a software environment for integrated models of biomolecular interaction networks, Genome Res 13 (2003), 2498–2504.

Szklarczyk

Morris

J.H.

Cook

Kuhn

Wyder

Simonovic

Santos

Doncheva

N.T.

Roth

Bork

Jensen

L.J.

and von Mering

, The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible, Nucleic Acids Research 45 (2017), D362–D368.

Koschützki

and Schreiber

, Centrality analysis methods for biological networks and their application to gene regulatory networks, Gene Regulation and Systems Biology 2 (2008), GRSB. S702.

Bader

G.D.

and Hogue

C.W.

, An automated method for finding molecular complexes in large protein interaction networks, BMC Bioinformatics 4 (2003), 2.

10.

Bindea

Mlecnik

Hackl

Charoentong

Tosolini

Kirilovsky

Fridman

W.H.

Pages

Trajanoski

and Galon

, ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks, Bioinformatics 25 (2009), 1091–1093.

11.

Hill

D.P.

Smith

McAndrews-Hill

M.S.

and Blake

J.A.

, Gene Ontology annotations: what they mean and where they come from, BMC Bioinformatics 9(Suppl 5) (2008), S2–S2.

12.

Hahn