Network-based analysis of differentially expressed genes in cerebrospinal fluid (CSF) and blood reveals new candidate genes for multiple sclerosis

Transcriptome data collection and processing

Gene expression profiles in both CSF cells and PBMCs were obtained from the ArrayExpress Database under the accession number of EMTAB- 69 based on the Human Genome 133 plus 2.0 arrays (Brynedal et al., 2010). Accordingly, this study consisted of 26 multiple sclerosis patients, of whom 12 and 14 patients were sampled during relapse and remission, respectively. The MS patients were selected from a large cohort of newly diagnosed MS patients, and none of the patients had ever received immunomodulatory drugs. Control population included 18 subjects with other neurological diseases to assess MS specific transcriptome. The microarray raw data were converted to gene expression values using the RMA algorithm by the affy package within R software (Gautier et al., 2004). After preprocessing, each expression profile containing 54, 675 probe sets that ones with less discriminative power were removed according to the measurement of overall variance by the varFilter function using the genefilter package from the Bioconductor project within R software (Gentleman et al., 2011). After the preprocessing, a total of 27,336 probe sets from each sample were used for further analysis. To identify differential expression of the selected probes, the limma package in R software was used to perform the moderated t-test (Smyth, 2005). Where a gene had more than one probe on the microarray, the average expression value of all the related probes was used to estimate expression level of the gene.

Interactome data

The human PPI network was gathered from four major IMEx (Orchard et al., 2007) public databases: IntAct (Kerrien et al., 2012), MINT (Ceol et al., 2010), DIP (Xenarios et al., 2002) and InnateDB (Lynn et al., 2008). Indeed, public PPI databases which only stored experimentally verified interactions used to eliminate possible spurious interactions and avoid misleading conclusions. Our recent study showed IMEx databases (especially IntAct and DIP databases) had a greater number of significant correlations for their proteins’ topological features than the all other paired comparisons between BIND, HPRD, MINT, IntAct and DIP databases (Safari-Alighiarloo, Taghizadeh & Rezaei-Tavirani, 2015).

QQPPI networks construction and topological analysis

To construct Query-Query PPI (QQPPI) networks, the differentially expressed genes in CSF (MS vs. control) and PBMCs (relapse vs. remission) were separately located on human PPI network. QQPPI networks only consistent of the query genes as the nodes and direct interactions among them. The subnetworks of QQPPI were constructed and visualized by Cytoscape software (Shannon et al., 2003). Centrality parameters of QQPPI networks were analyzed using the Cytoscape and CentiBin softwares (Junker, Koschützki & Schreiber, 2006). The following parameters were calculated to determine biologically significant nodes (Zhang, 2009). Degree: the number of links to a given node. The Betweenness centrality of node v is calculated as: (1)${\displaystyle C_B\left(v\right)={\sum }_{s\ne t\ne v\in V}\frac{\rho st\left(v\right)}{\rho st},}$ where the number of all shortest paths between node s and t regarded as ρst, and the number of shortest paths which passing through a node v out of ρst regarded as ρst(v). Indeed, this formula represents the ratio of the number of shortest paths passing through node v to the number of all shortest paths between s and t. The current flow betweenness centrality of a node v is the average of the current flow over all source–target pairs. Closeness centrality is defined as the reciprocal of the total distance from a node v, to all other nodes. Therefore, high values of closeness should indicate that all other nodes are in proximity to node v. (2)${\displaystyle C_C\left(v\right)=\frac{1}{\sum _{u\in v}\mathrm{dis}\left(u,v\right)}.}$ The centroid value is the most complex node centrality index and is computed by focusing the calculus on couples of nodes (v, w). The centroid value of an individual node ‘v’ is calculated by considering the number of nodes that have minimum shortest path which are closer to ‘v’ than ‘w’. A node v with the highest centroid value is the node with the highest number of neighbors separated by the shortest path to v, (3)${\displaystyle C_{cen}\left(v\right)=\min \left(f\left(v,w\right)\right),}$ where $f\left(v,w\right)=yv\left(w\right)-yw\left(v\right)$ and yv(w) denotes the number of nodes that are closer to v than w. Eigen vector centrality assign the relative significance of all nodes in the network by weighting connections to highly important nodes more than connections to nodes of low importance. (4)${\displaystyle {\lambda C}_{IV}={AC}_{IV},}$ where C_IV donates the Eigen vector and λ donates the Eigen value.

Hub and bottleneck nodes were extracted from the networks in two steps; (1) In the networks, nodes with degree greater than or equal to the sum of mean and twice the standard deviation (S.D.), i.e., mean + 2*S.D. of the degree distribution, were considered as hubs (Ray, Ruan & Zhang, 2008). (2) We defined bottlenecks as the proteins that were in the top 5% in terms of betweenness centrality.

Identification and annotation of functional modules

Network clustering was implemented by Clustering with overlap neighborhood expansion (ClusterONE) algorithm in order to identify the connected regions within the networks with possible overlap (Nepusz, Yu & Paccanaro, 2012). The modules were identified to have a minimum density of >0.05 and a degree of >5. A cluster with a p-value of <0.05 was determined to be a module. The functional meaning for identified modules was further explored, and they considered as candidate functional modules if their genes were significantly enriched in the biological process of Gene Ontology (GO) annotation or KEGG pathway.

Identification of complexes containing clique

CFinder software was applied to extract biologically meaningful protein complexes from the PPI networks (Adamcsek et al., 2006). CFinder (http://www.cfinder.org/) was downloaded and implemented locally. Cliques with 3 nodes and 4 nodes (3-cliques, 4-cliques) were identified in the QQPPI networks by this software. The cliques were searched against CORUM database (Ruepp et al., 2010) to find significant protein complexes. Then, all the proteins associated with a specific complex were identified using the in house algorithm. Complexes containing 3 or more query proteins, as a cut-off, were listed in this study.

Functional enrichment analysis

An enrichment analysis was performed using Functional Annotation Chart in DAVID bioinformatics. To determine functional modules, only the enriched GO terms and pathways with p-values < 0.05 were considered significant. Furthermore, Cytoscape Enrichment Map plugin was used to visualize significant terms enriched in entire networks by following parameters: p-value cut-off = 0.001, q-value cut-off = 0.05, overlap coefficient cut-off = 0.6 (Merico et al., 2010).

Expression analysis

We used the Limma package to analyze gene expression profile, E-MTAB-69, for comparison of four transcriptomes in MS (CSF: MS vs. controls and relapse vs. remission, PBMCs: MS vs. controls and relapse vs. remission). There were 3,062 genes with FDR < 0.05 whose expression was different in the CSF of MS patients as compared to the controls, but none in the respective PBMCs comparison. The number of up and down regulated genes was 1,080 and 1,982, respectively. On the contrary, when MS patients in relapse to those in remissions were compared, 1,163 differential expression genes with FDR ≤ 0.1 were seen in PBMCs, but none in the CSF. The number of up and down regulated genes was 301 and 762, respectively. The full lists of annotated differentially expressed probe sets are shown in Table S1 for the MS vs. control comparison in CSF cells, and in Table S2 for the relapse versus remission in PBMCs cells.

Networks’ topological analysis

We used only direct interactions of differentially expressed genes to construct QQPPI networks. The CSF PPI network consisted of 1,440 nodes and 3,500 edges and PBMCs PPI network involved 483 nodes and 941 edges. Topological features were processed to characterize the biology network from the random network. The power law of node degree distribution is one of most important criteria (Maslov & Sneppen, 2002; Zhu, Gerstein & Snyder, 2007). The distribution of node degree approximately followed power law distributions, where P(k) is a distribution of node degree, k is a degree and λ is a degree exponent, with λ = 1.94 and λ = 2.08 for CSF and PBMCs networks, respectively, and Fig. 1 indicates that the QQPPI networks were scale-free. The hubs and bottlenecks were extracted from the QQPPI networks by the criterion described in the method section (Table 1). Besides, we calculated four others centrality measurements involving closeness centrality, centroid value, Eigen vector centrality and current flow betweenness centrality and identified more central genes in the networks. The list of all nodes and their centrality measurements are prepared in Tables S3 and S4 for CSF and PBMCs in which candidate markers have been highlighted. The graphical structure of CSF and PBMCs PPI networks containing 5% top central genes are represented in Fig. 2.

Modularity analysis

The ClusterONE algorithm was selectively implemented on CSF and PBMCs networks to mine the functional modules, which may reveal a lot of hidden biological significant processes. We further performed GO and pathway analysis using DAVID tool to characterize functional modules; 14 and six functional modules were discovered for CSF and PBMCs (p-value < 0.05), respectively. In the case of CSF, enriched modules were relevant to the comparison of MS versus controls in which modules correlated remarkably with many immune-related pathways such as, cytokine–cytokine receptor interaction, chemokine signaling pathway, Toll-like receptor signaling pathway, T cell receptor signaling pathway and Hematopoietic cell lineage. Further to them, some modules were enriched for apoptosis, p53 signaling pathway, MAPK signaling pathway, Hedgehog signaling pathway and Fc gamma R-mediated phagocytosis. The other major enriched pathways in modules included focal adhesion, cell cycle, endocytosis, gap junction, tight junction, ECM-receptor interaction, regulation of actin cytoskeleton (Table 2).

For PBMCs, enriched modules corresponded to the comparison of relapse versus remission in which the majority of enriched pathways contributed to immune-related pathways like antigen processing and presentation, primary immunodeficiency, RIG-I-like receptor signaling pathway, Toll-like receptor signaling pathway and cytosolic DNA-sensing pathway. TGF-beta signaling pathway and p53 signaling pathway were the two noticeable signaling pathways in modules. The last enriched pathways were spliceosome, proteasome, ubiquitin mediated proteolysis and cell cycle (Table 3).

Identification of cliques and complexes

CFinder software was implemented to identify several 3-cliques and 4-cliques in the QQPPI networks. The corresponding complexes were retrieved from the CORUM database and shown in Table 4. In the case of CSF (MS vs. controls), these complexes mediated various biological functions such as protein processing (proteolytic), proteasomal degradation, stress response, protein binding, protease activator (ID: 32, 192 and 193), DNA conformation modification, transcription repression, protein modification by acetylation, deacetylation (ID:58), DNA conformation modification, transcription repression and posttranscriptional control (ID:105, 995,996 and 974), mitotic cell cycle and protein modification (ID:310 and 313), chromosome segregation/division (ID:1464), cell junction (ID:1839), actin cytoskeleton organization and biogenesis (ID:3008), ribosome biogenesis (ID:3055), protein modification and cellular signaling (MAPKKK cascade (ID:5909 and 5886).

For PBMCs (relapse vs. remission), the identified complexes involved in many biological processes like protein processing (proteolytic), proteasomal degradation, stress response, protein binding, protease activator (ID: 38, 39, 181, 191, 192, 193, 194), RNA processing and RNA binding (ID: 351 and ID: 1181), protein targeting, sorting and translocation, protein transport and homeostasis (ID:437), protein kinase (ID:5199), NIK-I-kappaB/NF-kappaB cascade and cytokine activity (ID: 5269), apoptosis (ID: 5473, 5860, 5859, 5799 and 5800).

Functional enrichment analysis of the networks

To gain a full view of the networks potential functions, the networks’ nodes were annotated using the Functional Annotation Chart in DAVID and visualized using the Enrichment Map plugin in Cytoscape. As shown in Fig. 3, each node represented one functional annotation term. Nodes with more enriched genes were larger. Edge width was indicated the extent of overlapping between these categories (overlap coefficient cut-off 0.6). In case of CSF (MS vs. control), the ten most enriched entries in Gene Ontology (GO) biological processes were GO:0009611∼Response to wounding, GO:0006955∼Immune response, GO:0007242∼Intracellular signaling, GO:0010604∼Positive regulation of macromolecule metabolic processes, GO:0007049∼Cell cycle, GO:0002684∼Positive regulation of immune system processes, GO:0042981∼Regulation of apoptosis, GO:0006954∼Inflammatory response, GO:0016044∼Membrane organization, GO:0008283∼Cell proliferation. In case of PBMCs (relapse vs. remission), ten most enriched terms included GO:0043068∼Positive regulation of programmed cell death, GO:0006974∼Response to DNA damage stimulus, GO:0043065∼Positive regulation of apoptosis, GO:0006917∼Induction of apoptosis, GO:0044265∼Cellular macromolecule metabolic processes, GO:0033554∼Cellular response to stress, GO:0022402∼Cell cycle processes, GO:0010033∼Response to organic substance, GO:0002684∼Positive regulation of immune system process, GO:0051249∼Regulation of lymphocyte activation.

Mining and identification disease markers in modules and complexes

We screened more central nodes in CSF and PBMCs to investigate detailed analysis about their association in functional modules and complexes. Interestingly it was found that eight genes (STK4, RB1, CDKN1A, CDK1, RAC1, ARRB2, ARRB1, FN1) were located in functional modules and 15 genes (RBBP4, ILF2, RPS16, PSMA3, EED, EZH2, CDK1, CDKN1A, CTNNB1, SDCBP, SRPK1, TUBA1A, YBX1, SFN, YWHAG) were associated in complexes in CSF-QQPPI network. In case of PBMCs-QQPPI network, four genes (CDK2, PSMA1, IKBKE and MYC) were located in functional modules and genes (PSMA1, PSMA2, PSMA7, PSMB3, PSMD3, HNRNPM, FAS, ACTB, CDC37, HSP90AB1, MAP3K3, IKBKE, CASP8) were associated in complexes. Besides, it was obvious from the mining of functional modules and complexes that more central genes incorporated in driving pathways of MS pathogenesis (Table 5, also see Tables 2, 3 and 4). The more central genes and their expression values in CSF and PBMCs networks are illustrated in Fig. 4 and candidate markers are represented in Fig. 5. Since these differentially expressed genes in microarray dataset corresponded to topologically significant nodes in PPI networks, which have functional importance because of their involvement in functional modules and complexes, they called as candidate disease markers in our study.