UBC Faculty Research and Publications

Meta-analysis of gene coexpression networks in the post-mortem prefrontal cortex of patients with schizophrenia… Mistry, Meeta; Gillis, Jesse; Pavlidis, Paul Sep 26, 2013

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-12868_2013_Article_3409.pdf [ 669.93kB ]
JSON: 52383-1.0223138.json
JSON-LD: 52383-1.0223138-ld.json
RDF/XML (Pretty): 52383-1.0223138-rdf.xml
RDF/JSON: 52383-1.0223138-rdf.json
Turtle: 52383-1.0223138-turtle.txt
N-Triples: 52383-1.0223138-rdf-ntriples.txt
Original Record: 52383-1.0223138-source.json
Full Text

Full Text

RESEARCH ARTICLE Open AccessMeta-analysis of gene coexpression networks inthe post-mortem prefrontal cortex of patientswith schizophrenia and unaffected controlsMeeta Mistry1, Jesse Gillis2 and Paul Pavlidis3,4*AbstractBackground: Gene expression profiling of the postmortem human brain is part of the effort to understand theneuropathological underpinnings of schizophrenia. Existing microarray studies have identified a large number ofgenes as candidates, but efforts to generate an integrated view of molecular and cellular changes underlying theillness are few. Here, we have applied a novel approach to combining coexpression data across seven postmortemhuman brain studies of schizophrenia.Results: We generated separate coexpression networks for the control and schizophrenia prefrontal cortex andfound that differences in global network properties were small. We analyzed gene coexpression relationships ofpreviously identified differentially expressed ‘schizophrenia genes’. Evaluation of network properties revealeddifferences for the up- and down-regulated ‘schizophrenia genes’ , with clustering coefficient displaying particularlyinteresting trends. We identified modules of coexpressed genes in each network and characterized them accordingto disease association and cell type specificity. Functional enrichment analysis of modules in each network revealedthat genes with altered expression in schizophrenia associate with modules representing biological processes suchas oxidative phosphorylation, myelination, synaptic transmission and immune function. Although a immune-functionenriched module was found in both networks, many of the genes in the modules were different. Specifically, adecrease in clustering of immune activation genes in the schizophrenia network was coupled with the loss of variousastrocyte marker genes and the schizophrenia candidate genes.Conclusion: Our novel network-based approach for evaluating gene coexpression provides results that converge withexisting evidence from genetic and genomic studies to support an immunological link to the pathophysiology ofschizophrenia.Keywords: Schizophrenia, Microarray, Gene coexpression network, Postmortem brainBackgroundSchizophrenia is a severe psychiatric disorder with anelusive etiology. Gene expression profiling of the post-mortem human brain has been frequently used as ameans to investigate patterns of molecular disruption inthe brains of patients with schizophrenia. One of themost common types of analysis applied to expressionprofiling data is differential expression; which is used toidentify over- or under-expressed genes associated with theillness. Candidate genes identified from expression profilingstudies in schizophrenia have implicated alterations indifferent cellular systems, including myelination, synaptictransmission, metabolism, and ubiquitination [1]. Thesefindings are not always replicated across studies, nor havethey been successfully integrated into a comprehensivebiological framework.In our previous work, we used a large combinedcohort to identify a meta-signature of genes which areconsistently differentially expressed in the prefrontal cor-tex of patients with schizophrenia [2]. The functionsreflected in these genes are diverse and the interactionsamong them are largely unexplored. Because gene* Correspondence: paul@chibi.ubc.ca3Department of Psychiatry, University of British Columbia, Vancouver, BC,Canada4Centre for High-throughput Biology, 177 Michael Smith Laboratories, 2185East Mall, University of British Columbia, Vancouver, BC, CanadaFull list of author information is available at the end of the article© 2013 Mistry et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the CreativeCommons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, andreproduction in any medium, provided the original work is properly cited.Mistry et al. BMC Neuroscience 2013, 14:105http://www.biomedcentral.com/1471-2202/14/105function is partly defined by interactions with othergenes (at the biochemical, physical interaction, geneticor regulatory levels), it is attractive to apply genecoexpression network analysis to aid in interpretation.In general, gene networks can be analyzed to identifyhigher-level features of gene-gene relationships based ongraph theoretic considerations such as node degree orclustering coefficient [3-5]. Evaluating the broader net-work structure allows us to detect modularity in thegraph, or groups of densely connected nodes with sparseconnections between groups. Characterization of these‘modules’ can convey useful information as they may beassociated with specific molecular complexes or func-tions, yielding hypotheses that would be difficult to as-certain based on a gene-by-gene analysis. It is importantto note that the terminology “gene coexpression net-work” refers to a sparse representation of the correlationstructure among genes, and that such networks are notamenable to straightforward interpretation in the waythat protein interaction or metabolic networks are.However, a key advantage of coexpression is that there isrelatively abundant data, so “condition-specific networks”can be constructed. Thus one can evaluate differences be-tween condition-specific networks to help elucidate systemslevel molecular dysfunction.Such coexpression network analyses have recentlybeen applied to a number of postmortem human brainexpression profiling datasets for examining general tran-scriptome patterns of the CNS [6], and to interrogatethe molecular basis of neuropsychiatric disease [7-10].Torkamani and colleagues [8] conducted a network ana-lysis by combining two independent schizophrenia ex-pression profiling datasets. Expression data was mergedacross control and schizophrenia cohorts and modulesof coexpressed genes were characterized according todisease characterization, cell type specificity and the ef-fects of aging. A more recent cross-cortical networkstudy was carried out by Roussos et al. [10] using controland schizophrenia samples across four different brain re-gions. Discrete modules of coexpressed genes displayedhigh preservation between control and schizophrenianetworks for all but one module. Brain regional differ-ences were assessed with an analysis of variance com-parison of module eigengene expression, with changesonly observed in the control network. Chen et al. [11]also explored networks using combined data fromschizophrenics and controls. Two modules were associ-ated with genes differentially expressed with diseaseacross the datasets; one which was specific to cerebellarcortex and the other identified across all brain regions.They did not report any differences in networks betweenschizophrenics and controls. Although Chen et al. usedfour data sets, they were not independent as three of thedatasets used samples from the same brain collection.In this study, we applied coexpression network analysisto seven independent gene expression studies of theprefrontal cortex to demonstrate, in agreement withprevious studies, an overall similarity in transcriptomeorganization between healthy controls and individualswith schizophrenia. We then examined network proper-ties of genes we previously reported to be differentiallyexpressed in schizophrenia [2] within each network toreveal features of these genes that are not observed withother functional gene groups or other brain-relateddisease gene sets. Finally, using a network clustering ap-proach, we found evidence for functional dysregulationof immune-related processes in schizophrenia. Our re-sults complement previous gene expression and geneticsevidence supporting an immunological aspect of thedisorder.ResultsWe constructed two gene coexpression networks; onerepresenting the control human prefrontal cortex andthe other representing the prefrontal cortex in schizo-phrenia (referred to as CTL and SZ, respectively). Theschizophrenia and control groups had no significant dif-ferences in age and PMI (Table 1). Sex differences wereassessed by use of a chi-squared test for equality ofproportion, and we observed no significant difference(p = 0.1). Brain pH was significantly different (t-test;p = 0.001). Each network was comprised of 12,582 genes(nodes), and 392,606 coexpression ‘links’ among them.The two networks had similar values in the average clus-tering coefficient (p > 0.1), but average shortest pathlength across nodes differed slightly (p < 0.01). Thesenetwork properties are summarized in Table 2. Both net-works exhibited a ‘heavy-tailed’ node degree distribution,with most of the genes interacting with few partners anda small proportion of genes displaying ‘hub’-like behav-iour interacting with many genes. In the literature, suchdistributions are sometimes described as ‘scale free’. Weused a linear regression of the log-scale node degree dis-tribution to examine this in our networks (Figure 1).Table 1 Summary of demographic variables acrosscombined cohortControl Schizophrenia P-valueNumber of Subjects 153 153Age 56.25 ± 20 55.27 ± 19 p = 0.67Sex 101 M : 52 F 113 M : 40 F p = 0.1Brain pH 6.5 ± 0.28 6.39 ± 0.29 p = 0.001PMI 21.95 ± 15.3 22.65 ± 15.2 p = 0.69F, female; M, male; PMI, post-mortem interval. There were 319 samplescollected across seven datasets of which 306 passed quality control analysis.The summary demographics (mean ± standard deviation) and t- test p-valuesfor group differences are shown for those subjects used in the analysis. Forsex we report the p-value generated from a chi-squared test for equalityof proportions.Mistry et al. BMC Neuroscience 2013, 14:105 Page 2 of 16http://www.biomedcentral.com/1471-2202/14/105While a linear fit explains over 80% of the variance innode degree distribution, (CTL R2 = 0.857; SZ R2 =0.872), the fit is not good at the extremes. Based onmore stringent criteria (which we endorse), our net-works are not ‘scale-free’ [12]. However, the ‘heavy-tailed’ nature of the node degree distributions in ournetworks is typical of other ‘biologically relevant’ net-works cited in the literature [13,14].The small differences in global network properties ob-served between CTL and SZ suggests that there is anoverall common coexpression structure of the prefrontalcortex. Fifty-seven percent of the edges (224,384 links)are the same in the two networks, much higher thanexpected by chance. The remaining 168,222 edges arenot shared between the two networks (Figure 2). Subtledifferences between the networks are also indicated by ahigher maximum node degree in SZ (935) than CTL(737), and the increased number of non-connectednodes in CTL (2356) compared to SZ (2288). These dif-ferences could indicate subtle biological differences be-tween the two networks, but are presumably at leastpartly due to the effect of noise.In addition to comparing average network propertiesacross the SZ and CTL networks to each other, we com-pared each separately to a node degree-matched randomnetwork (see Methods). For features based on connectiv-ity (i.e. shortest path length and clustering coefficient),we found the observed distributions of both networks tobe higher than compared to random networks. Shortestpath length displayed slightly higher values than foundin randomized networks (Figure 3A, C). Additionally,genes showed an increased clustering into local commu-nities compared to genes from a randomized networkwith identical degree distribution (Figure 3B, D). Thuswhile, the SZ and CTL networks are similar, they arealso clearly distinct from random networks with thesame node degree distribution.Table 2 Whole network properties of the control andschizophrenia brain networksControl SchizophreniaNon-connected nodes 2356 2288Maximum node degree 747 935Mean node degree 77 76Shortest path length 3.34 3.32Cluster coefficient 0.29 0.29log-log fit (R2) 0.857 0.872Number of modules 25 25Figure 1 Connectivity distribution of control and schizophrenia networks. The control brain network (A) and the schizophrenia brainnetwork (B) connectivity distribution on a log10-log10 scale. Plotted on the x-axis is the number of links versus the number of genes that havethe corresponding number of links on the y-axis.Figure 2 Comparing node degree between networks. To assessthe node degree differences between networks, values of the nodedegree in the control network were plotted against the number ofedges retained in the schizophrenia network. Each gene isrepresented by a dot in the plot, and all 12,582 genes were plotted.The presence of data points that deviate from the identity lineindicate differences in gene-to-gene connections between the twonetworks. A histogram is also provided (inset) to illustrate that thedistribution mean of the overlap is about fifty percent.Mistry et al. BMC Neuroscience 2013, 14:105 Page 3 of 16http://www.biomedcentral.com/1471-2202/14/105We next investigated network properties at the level ofgene groups, focusing on our previously identified meta-signature of genes differentially expressed in schizophre-nia [2]. The meta-signature of 94 genes (25 up-regulatedand 69 down-regulated) will be referred to as SZUPand SZDOWN, respectively. Network properties wereassessed for each gene set individually by taking an aver-age across all genes in the group. These results aresummarized in Table 3. For each gene set we computedthe average values for shortest path length, cluster coef-ficient and node degree and evaluated differences ob-served between the control and schizophrenia networks.In general, both gene sets had a low mean node degreewith respect to the network degree distribution of CTLand SZ, tending not to be ‘hubs’. For the SZUP gene set,we found higher node degree, shorter path length andan increased clustering coefficient in the SZ network,though these differences were not statistically significant(Wilcoxon- Rank Sum test p > 0.05). Conversely, theSZDOWN gene set exhibited a decreased node degree,larger path length and a lower clustering coefficient(Wilcoxon- Rank Sum test p = 0.05) in the SZ network.Thus, the properties of each gene set displayed an appar-ent trend between networks and those trends were ob-served to be opposite between the two gene sets.The trends observed for SZUP and SZDOWN be-tween the two networks were small and only marginallysignificant. To evaluate whether each of the individualgene set values were unusual in the networks weimplemented three different methods of control. A firstcontrol was supplied by comparing observed networkmeasures for SZUP and SZDOWN to a backgrounddistribution of 1000 randomly selected gene setsof matched size and node degree (see Methods). Thedifference between the observed values and backgroundwas assessed by computing z-scores and p-values, asreported in Table 3. Our strongest result was for theclustering coefficient of both gene sets. The p-values forthe SZUP clustering coefficient indicate that the highvalue in SZ is significant when compared to a back-ground distribution (p = 0.005). For SZDOWN, the highclustering coefficient in CTL showed a trend differenceFigure 3 Comparison to random network distributions. For each network, we generated a corresponding random network by swappingedges and maintaining the same node degree distribution. (A, C) Shortest path length distribution of real networks are shifted slightly higherthan corresponding random network distributions, but distributions between CTL and SZ are similar. Grey histograms reflect values from therandom network, and black histograms represent the real network data. (B, D) Genes cluster into local communities with high number ofinterconnections compared to corresponding random networks. Black dots represent real network data and grey dots represent randomnetwork data.Mistry et al. BMC Neuroscience 2013, 14:105 Page 4 of 16http://www.biomedcentral.com/1471-2202/14/105compared to the background distribution (p = 0.06). To-gether, these results converge to highlight that the geneneighbouring meta-signature gene sets are unusuallyhighly coexpressed, with the SZUP genes displaying thisproperty in the SZ network and the SZDOWN genestending to display this in the CTL network.A second control was applied to examine whether ornot the properties observed for SZUP and SZDOWNare a feature of other functionally grouped sets of genes.This is a more stringent control than simply comparingto random gene sets, because we are interested in prop-erties of our genes that are unusual in schizophreniacompared to other functionally related groups of genes.We created 3,230 different functional gene sets based onGO terms and their associated genes. We assessed net-work measures for each functional gene set and com-pared the resulting z-scores to values for SZUP andSZDOWN (Figure 4). For the clustering coefficient, thez-score for SZUP is more distinguishable from GOgroup values in the SZ network compared to CTL. Theopposite is true for the SZDOWN z-score values. Thus,the clustering coefficient is a unique property of ourmeta-signature genes not observed with other functionalgene sets.We next evaluated whether our meta-signature genesets share network properties with gene sets associatedwith other brain-related disorders. We assembled genesets for five different disorders, mostly based on findingsfrom genetic association studies. For each disease geneset, z-scores were computed based on a background dis-tribution and compared against the GO group z-scoresin each network (Additional file 1: Table S1). Ofparticular interest are the results observed for clusteringcoefficient in the two networks. Interestingly, theAlzheimer’s disease gene group (Figure 4, red arrows)exhibited strikingly similar properties to SZUP in bothnetworks despite having only one overlapping gene. Not-ably, the Parkinson’s disease gene group (Figure 4, bluearrows) follows a similar but more subtle trend asSZDOWN.One concern is that these features are particular prop-erties of the data sets we analyzed. In the absence of suf-ficient independent data, we assessed the robustness ofthe network measures observed for SZUP andSZDOWN using a jackknife procedure. In this process,we removed one of the seven datasets and regeneratedaggregate CTL and SZ networks on the remaining six,for each study in turn. This yields seven pairs of jack-knife networks. For each jackknifed network, the averageshortest path length and clustering coefficient was com-puted for SZUP and SZDOWN and values were com-pared between networks (Additional file 2: Figure S1).For SZUP, we observed a general agreement of increas-ing clustering coefficient and consistently decreasingpath length between CTL and SZ across all iterations.For SZDOWN, we found that only the clustering coeffi-cient effects were robust to removing single data sets;the path length results proved to be more sensitive.Taken as a whole, these results are consistent with subtlenetwork property differences for the SZUP andSZDOWN genes between the two networks.Our analysis to this point examined either the entirenetworks or used supervised approaches to select sets ofgenes for analysis. We complemented this with an un-supervised method based on clustering [5]. This analysiswas motivated in part by the observation that the meta-signature gene sets showed significant modularity differ-ences between CTL and SZ. We hypothesized that theremight be additional differences, beyond the parts of thenetwork involving the meta-signature genes, or that thisanalysis might uncover additional features of the meta-signatures. Clustering resulted in 25 modules of varyingsizes, in each network. An overlap comparison betweenthe two sets of modules revealed strongly “matching”modules for 15 modules (p < 0.001, hypergeometric test;Figure 5), prompting further characterization of thesemodules according to disease association, cell-type spe-cificity and functional roles.We identified five modules in each network whichdisplayed the most significant association with genesdifferentially expressed in schizophrenia (Wilcoxonrank-sum test; p < << 0.001), summarized in Table 4(a list of all network modules is provided in Additionalfile 3: Table S2). In order to functionally characterizecoexpression modules we cross-referenced the modulegene lists with published lists of genes encoding markersTable 3 Schizophrenia gene set network propertiesUp-regulated(25)Down-regulated(69)CTL SZ CTL SZNode degreeMean 63.9 83.5 127.4 106.2Non-interacting nodes 2 2 7 3Edges (within gene set) 12 23 129 144Shortest PathMean 3.28 3.14 3.31 3.48Random gene set comparisonZ-score −0.58 −0.91 3.15 4.46p-value 0.23 0.15 0 0Cluster CoefficientMean 0.35 0.38 0.32* 0.27*Random gene set comparisonZ-score 1.11 2.47 1.51 −0.607p-value 0.14 0.005 0.06 0.28*Difference is significant between CTL and SZ at p = 0.05.Mistry et al. BMC Neuroscience 2013, 14:105 Page 5 of 16http://www.biomedcentral.com/1471-2202/14/105of different CNS cell types (Table 5; Additional file 4:Table S3) and examined functional roles by GO enrich-ment (Additional file 5: Table S4 and Additional file 6:Table S5). The top schizophrenia-associated module ineach network showed a significant enrichment of neuronalmarkers as well as genes belonging to GO categoriesrelated to the electron transport chain and oxidativephosphorylation (Module CTL24 and Module SZ3). A largeproportion of the SZDOWN genes were also containedin this module, in both networks. Both networks alsoidentified a module enriched for genes belonging toimmune-related GO categories (Module CTL1 andModule SZ15). These modules were highly similar betweennetworks, with a 71 percent overlap in gene membershipand both associated with schizophrenia disease status. Celltype enrichment of this module in the control networkidentified 44 astrocyte marker genes (at low stringencythreshold). Interestingly, the corresponding module in theschizophrenia network was reduced to eleven astrocytemarker genes. Also, we observed three SZUP genesFigure 4 Comparison of gene set properties to functional GO groups. Histograms represent z-score distributions for cluster coefficient (A, B)and shortest path length (C, D) computed across 3,230 different GO groups in the control and schizophrenia networks. Z-scores represent thedifference between the mean value of the network measure of the GO group compared to the mean of random gene sets of the same size andmatched node degree. Dashed lines plotted represent the z-score obtained for the down-regulated meta-signature gene set and solid linesrepresent the z-score obtained for the up-regulated gene set. Z-score values for selected brain-related disease gene sets are displayed ascoloured arrows on the x-axis. ALZ Alzheimer’s disease; PD Parkinson’s disease.Mistry et al. BMC Neuroscience 2013, 14:105 Page 6 of 16http://www.biomedcentral.com/1471-2202/14/105(BAZ1A, TMEM176A, and FTL) in this module in thecontrol network, but not in the schizophrenia network.The next two sets of modules that showed significantassociation with schizophrenia were also notable in highenrichment of cell-type marker genes. Module CTL25and Module SZ20 were both significantly enriched foroligodendrocyte markers genes (e.g. OLIG2, MBP,MAG) with a 61 percent overlap in total gene member-ship, and each contained two SZUP genes. Highlyranked GO terms from functional enrichment analysisof each module included the presence of terms such as“ensheathment of neurons” and “regulation of action po-tential” suggesting this is was a myelination-relatedmodule. In contrast, modules CTL20 and SZ9 (whichhave 71 percent overlap in gene membership) were sig-nificantly enriched for astrocyte marker genes (Table 5;p < << 0.001). The common functional role suggested byGO enrichment centered around glutamine metabolism(e.g. GAD1, GAD2; see Additional file 5: Table S4 andAdditional file 6: Table S5). SZUP genes also identifiedwith this module, with a slightly higher number of genesidentified in the SZ network (Table 3; p < < 0.001).The final top disease-associated modules from eachnetwork, CTL18 and SZ2, were each associated with dif-ferent cellular processes. Module CTL18 was enrichedfor neuronal marker genes and functional enrichmentanalysis revealed an association with genes related toneurotransmitter secretion. Furthermore, ten of theSZDOWN genes overlapped with this module. In theschizophrenia network, twelve SZDOWN genes (none ofwhich overlap with the previous ten) were identifiedin Module SZ2. Module SZ2 is also neuron markerenriched but is associated with genes involved inubiquitination. The two modules showed very little over-lap with each other in terms of gene membership (withonly six genes in common) and each highlighted differ-ent functional roles, yet both modules exhibited an asso-ciation with schizophrenia in their respective networks.We next evaluated the impact of covariates on our re-sults by using expression changes known to be associ-ated with age, brain pH and sex. Gene lists for thesefactors were compiled from a previous study of healthycontrol postmortem brain [16], and hypergeometricprobabilities were computed to evaluate the significanceof overlap with each module (Additional file 7: Table S6).No significant overlap was observed with the sex genes ineither network. The age and pH genes displayed enrich-ment across the top modules which was mostly consistentFigure 5 Heatmap comparison of modules between networks.Each network was clustered using WGCNA-based methods. Moduleswere compared between networks by computing the number ofoverlapping genes and a percent overlap was computed by takingthe module of smaller size as the denominator. For each modulepair the percent overlap is plotted with the scale provided.Table 4 Disease characterization of coexpression modulesin each networkControl NetworkModule Theme ModuleNo.ModuleSizeSZUPOverlapSZDOWNOverlap(25 genes) (73 genes)OxidativephosphorylationCTL24 1996 0 34***GlutaminemetabolismCTL20 504 7*** 0Myelination CTL25 749 2 0Immune response CTL1 329 3** 0SynaptictransmissionCTL18 291 0 10***Schizophrenia NetworkModule Theme ModuleNo.ModuleSizeSZUPOverlapSZDOWNOverlap(25 genes) (73 genes)OxidativephosphorylationSZ3 1144 0 34***GlutaminemetabolismSZ9 717 12*** 1Myelination SZ20 711 2 0Immune response SZ15 198 0 0Ubiquitination SZ2 1424 0 12The top five modules in each network were characterized by enrichment ofgenes that are differentially expressed in schizophrenia. Module themes wereassigned based on Gene Ontology enrichment analysis. Disease effect p-valueswere computed by entering the t-statistic for the disease effect of each geneinto a Wilcoxon rank-sum test by module. For all modules listed the diseaseeffect p-values are p < << 0.001 and values can be found in the supplement.SZUP and SZDOWN are the up- and down-regulated schizophrenia gene setspreviously identified in our meta-analysis [2]. **p < 0.01; ***p < < 0.001.Mistry et al. BMC Neuroscience 2013, 14:105 Page 7 of 16http://www.biomedcentral.com/1471-2202/14/105between the control and schizophrenia networks. How-ever, differential enrichment was observed for the immuneresponse module, a trend similar to what we observedwith SZUP genes and astrocyte marker genes. Expressiondata for the age and pH genes from the immune responsemodule were plotted and fold change values were com-puted to examine the extent to which these genes mightbe affecting the network changes we observed (Additionalfile 8: Figure S2). We found expression to be variablewithin each cohort and differences in mean expression be-tween cohorts were very small and not significant. Thus,it seems unlikely that there are large confounding effectsof age and pH driving the changes we see in coexpressionclustering between the two networks. We also addressedmedication effects by cross-referencing lists of gene ex-pression changes associated with lifetime antipsychoticuse, to the top modules identified in our network analysis.The SMRI Online Genomics Database (https://www.stanleygenomics.org/) provides gene lists for several demo-graphic variables which have been independently assessedusing only the diseased subsets of each dataset. Fromthe database we extracted significant gene lists (p < 0.001;FC > 1.2) pertaining to the effects of lifetime antipsy-chotics (69 genes) in subjects with schizophrenia. In thecontrol network, a total of 34 medication-associatedgenes overlapped with our top modules, with the lar-gest number found in the immune response module(Additional file 7: Table S6). In the schizophrenia network,roughly the same number of genes are re-distributedacross the modules, but a significant overlap with theimmune response module remains.DiscussionOur network-based approach for evaluating genecoexpression provides a novel assessment of coexpressionpatterns across seven large schizophrenia microarraydatasets. We implemented a rank aggregation approachfor network analysis revealing interesting patterns ofmolecular connectivity in the control and schizophreniapostmortem human brain. Overall, the two coexpressionnetworks were very similar to one another. This is consist-ent with existing findings from network analysis inschizophrenia [8,10,11]. The control and schizophrenianetworks shared a similar node degree distribution, andaverage values of path length and clustering coefficienttaken across all nodes in the network were not signifi-cantly different. However, closer inspection revealed dif-ferences of potential biological significance.To evaluate differences in gene-gene connectivity be-tween networks, we initially focused on the networkproperties of 95 differentially expressed ‘schizophreniagenes’ as reported in our previous study of these samedata sets [2]. This gene list was divided into two groups:1) genes which are up-regulated in schizophrenia and 2)genes which are down-regulated in schizophrenia. Weexamined the network properties of each gene set withinthe control and schizophrenia networks and identifieddistinguishing features of our ‘schizophrenia genes’. Theclustering coefficient, a measure which gives us insightinto the community structure of nodes, proved to be aninteresting characteristic of both gene sets. Importantly,we applied control protocols to demonstrate that thisdifferential coexpression among the neighbourhood of‘schizophrenia genes’ is not observed with other groupsof “functionally related” genes and most other brain-related disease gene groups.We also performed an assessment of modularity acrossall nodes in each network. Our results were comparableto two module coexpression analyses previouslyconducted in schizophrenia. The top disease-associatedmodule in both networks was enriched for genes in-volved in oxidative phosphorylation and energy produc-tion This is consistent with results from Torkamani andcolleagues [8] in which a combined network was gener-ated from two expression datasets, both of which wereTable 5 Cell-type marker enrichment of coexpressionmodules in each networkControl NetworkLow stringency threshold ( > 4-fold)Module Theme ModuleNo.Oligodendrocyte Neuron AstrocyteOxidativephosphorylationCTL24 43 253*** 48GlutaminemetabolismCTL20 23* 9 166***Myelination CTL25 115*** 16 42Immune-related CTL1 7 10 44***SynaptictransmissionCTL18 7 85*** 9Schizophrenia NetworkLow stringency threshold ( > 4-fold)Module Theme ModuleNo.Oligodendrocyte Neuron AstrocyteOxidativephosphorylationSZ3 23 181*** 22GlutaminemetabolismSZ9 41*** 20 203***Myelination SZ20 109*** 16 28Immune-related SZ15 3 4 11Ubiquitination SZ2 34 157*** 36For the top disease-associated clusters in each network we report the numberof genes that overlap with published lists of cell-type marker genes foroligodendrocytes , neurons, and astrocyte marker genes provided by [15]. Cell-type enrichment for all modules in each network and results for highstringency lists (> 10-fold) can be found in Additional file 4: Table 1.Hypergeometric probabilities were computed to evaluate significance ofoverlap. *p < 0.05; ***p < < 0.001.Mistry et al. BMC Neuroscience 2013, 14:105 Page 8 of 16http://www.biomedcentral.com/1471-2202/14/105used in our study. Moreover, our finding fits within abody of evidence suggesting mitochondrial dysfunctionand defects in brain metabolism leading to oxidativestress in schizophrenia [17,18]. An oligodendrocytemarker enriched module was also identified in bothCTL and SZ networks of our study and supported withthe association of relevant GO terms related to myelin-ation. In agreement, Roussos et al. [10] identified theirtop schizophrenia-associated module as highly enrichedfor oligodendrocyte marker genes as well as for genes asso-ciated with cytoskeleton rearrangement, axonal guidanceand synaptogenesis. We note that the data used by Roussoset al. includes some samples contained in the Haroutuniandataset we used, albeit preprocessed in a different way(Table 6). The association of a myelination-related modulewith schizophrenia is in line with a wide range of reportedwhite matter abnormalities linked to the illness [23] andgenetic studies that have contributed a number of myelinand oligodendrocyte –related genes as candidate genes(e.g. APOD, PLP1, MAG) [24-26]. Torkamani and col-leagues also identified an oligodendrocyte/myelin-relatedmodule in their network, but they did not observe any asso-ciation with genes differentially expressed in schizophrenia.In our analysis, an ‘immune’ module consistentlyappeared in both networks. While many genes in thismodule were conserved between the two networks, wefound a number of differences suggesting alteration ofimmune-related processes in schizophrenia. In the con-trol network, the immune response module was muchlarger in size and contained four times as many astrocytemarker genes (at low-stringency threshold) than in theschizophrenia network. A list of microglia marker geneswere obtained from Bedard et al. (2007) [27] and alsocross-referenced with top modules from each network.The numbers of overlapping genes with each modulewere few and did not change between the two networks.Although microglia are considered the resident macro-phages of the brain providing the main arm of immunedefense in the CNS, much evidence suggests that astro-cytes also play an important role in the local regulationof immune reactivity [28]. In the control network im-mune module, 31 of the 44 overlapping astrocyte markergenes were coexpressed with one or more immune acti-vation genes (Figure 6). Moreover, most of those astro-cyte marker genes do not appear in the schizophrenianetwork module as a result of the loss of the immuneactivation genes. Torkamani et al. also identified an im-mune module, however, similar to the myelination mod-ule they failed to find any association with differentiallyexpressed schizophrenia genes. Importantly, in our ana-lysis three of the SZUP genes, specifically FTL, BAZ1Aand TMEM176A, were identified in this module in CTLbut not in SZ.BAZ1A encodes a subunit of the chromatin assemblyfactor (ACF), which together with other proteins com-prises the chromatin remodeling complex. Although it isnot directly associated with the immune system, wefound that it was coexpressed with a number of immuneactivation genes (Figure 6), suggesting a possible role inthe transcriptional regulation of immune related genes.TMEM176A was also coexpressed with immune-relatedgenes, complement component 4A (C4A) and complementcomponent 4B (C4B), both of which are also not found inthe SZ immune module (Figure 6). TMEM176A encodes atransmembrane protein, which together with TMEM176Bwhen overexpressed has been shown to block dendritic cellmaturation in rats [29]. Dendritic cells (DC) are immunecells found in most major organs, and their ability to regu-late immunity is dependent on DC maturation. The brainhas long been considered devoid of DC in the absence ofinflammation, with microglia charged with many functionalattributes commonly ascribed to DC. However, recent evi-dence has illustrated that DC are found in various tissuereservoirs within the steady-state CNS and are also poten-tial players in brain immune surveillance [30].Genes differentially expressed between schizophreniasubjects and healthy controls have also been identifiedthrough a recent combined analysis conducted on six ofthe same datasets used in our study [31]. Perez-Santiagoet al. identified 117 up-regulated and 43 down-regulatedTable 6 Datasets used in coexpression network analysisDataset Reference MicroarrayPlatformBrain region(s) No. of SubjectsCTL:SZStanley Bahn SMRI database HG-U133A Frontal BA46 31: 34Stanley AltarC SMRI database HG-U133A Frontal BA46/10 11: 9*Mclean HBTRC HG-U133A Prefrontal cortex (BA9) 26: 19Mirnics Garbett K. et al., 2008 [19] HG-U133A/B Prefrontal cortex (BA46) 6: 9*Haroutunian Katsel P. et al., 2005 [20] HG-U133A/B Frontal (BA10/46) 29: 31GSE17612 Maycox P. et al., 2009 [21] HG-U133 Plus 2.0 Anterior prefrontal cortex (BA10) 21: 26*GSE21138 Narayan S. et al., 2008 [22] HG-U133 Plus 2.0 Frontal (BA46) 29: 25SMRI, Stanley Medical Research Institute; HBTRC, Harvard Brain Tissue Resource Centre (Mclean66 collection).*Samples from these datasets have been used in previous coexpression network studies of schizophrenia [8,10].Mistry et al. BMC Neuroscience 2013, 14:105 Page 9 of 16http://www.biomedcentral.com/1471-2202/14/105probes of which 8 and 10 probes overlap with our up-and down- meta-signatures, respectively. The genesfrom Perez-Santiago et al. displayed a similar enrichmentpattern as the Mistry et al. signature genes in the topmodules in each of our networks. Specifically, the up-regulated genes showed a trend that coincides with trendsobserved with the astrocyte marker genes and SZUP genes.In the control network, the immune module contains 44up-regulated genes and in the schizophrenia network thenumber of genes in the immune module is reduced to 19genes. These findings provide additional support to our dis-cussion on genes up-regulated in schizophrenia and alter-ation of immune-related processes.The immunological link to the pathophysiology ofschizophrenia suggested from our network analysis isnot a new concept. Linkage and GWAS support an asso-ciation of a broad section of markers in the major histo-compatibility complex (MHC) region at 6p21.33 withschizophrenia [32,33] and many immune-related genes,genetic variants and haplotypes are also implicated inschizophrenia [34-36]. Several studies of gene expressionin the postmortem PFC have reported alterations inimmune and stress-response genes in subjects withschizophrenia [37-39]. Additionally, the investigation ofgene expression in peripheral blood mononuclear cells(PBMCs) of drug-naïve schizophrenia subjects also sup-ports alteration of immune-related processes [40-42].Using our panel of schizophrenia genes found frommeta-analysis we were able to identify unique features ofthe coexpression network and highlight relevant areas ofdysfunction which may contribute to the pathophysi-ology of schizophrenia. However, our interpretation ofthe network model is based on GO enrichment and fur-ther investigation at the individual gene level will pro-vide an explanation of higher resolution. We also soughtevidence to support the plausibility of schizophreniaFigure 6 Gene membership differences for the ‘immune’ module between networks. A) Visualization of the ‘immune’ module in thecontrol network is depicted by taking 28 immune activation genes (GO:002253) form the module and their associated connections within themodule. This resulted in a sub-module of 229 nodes and 672 edges. Immune activation genes are represented by the larger nodes on the outerring, astrocyte marker genes are represented by the smaller nodes on the inner ring and ‘schizophrenia genes’ are represented by diamondshaped nodes. Lighter colored nodes indicate the gene is not present in the schizophrenia network immune module while darker colored nodesindicate the gene is retained. B) and C) Bar graphs were used to illustrate the difference in cell-type marker membership for the ‘immune’module in the two networks (using the low-stringency list).Mistry et al. BMC Neuroscience 2013, 14:105 Page 10 of 16http://www.biomedcentral.com/1471-2202/14/105genetic risk variants driving the changes we observed.From the SZGene (http://www.szgene.org/) database wecompiled the top 45 most reliable genetic associationsbased on findings from a meta-analysis [43] . This genelist was cross-referenced with the top modules fromeach network and we found no overlap with the immunemodule in either network.One conclusion of our study is the effects on genecoexpression patterns in schizophrenia, like those on ex-pression levels, are subtle in the face of other sources ofvariability. This suggests that while seven datasets is amuch larger cohort than used in previous studies,our study would benefit from having additional data.Aggregating across a larger number of datasets has beenshown to result in networks more comparable to PPInetworks [44], and is likely to give more reliable interac-tions, but saturation of some properties was not reportedto be achieved until >20 networks were combined. To helpameliorate concerns of robustness, we used a jackknife ap-proach. The results from the jackknife analysis also addressconcerns regarding the reuse of datasets for both differen-tial expression [2] and coexpression. We established thatour findings are not overly sensitive to influence of anysingle data set. Also, for our null model comparisons weapplied stringent controls whereby node degree was con-trolled for. Because node degree is one of the most import-ant features of a network, it can drive numerous effects onother topological properties such as clustering coefficientand shortest path length. Thus controlling for node degreein generating null distributions is critical for avoiding falsepositives. Similarly, because we wished to identify networkfeatures which are associated with schizophrenia, we con-trolled for “generic” network features by comparing ourschizophrenia gene sets to the properties of other “func-tionally coherent” gene groups, not just randomly selectedgenes. This approach was motivated by recent work fromour group showing that generic multifunctionality effectscan strongly skew the interpretation of gene networks [45].As is the case with most postmortem brain studies inschizophrenia, our results should be interpreted in thecontext of several caveats. Samples used in this studywere taken from patients having lived with schizophre-nia for various lengths of time, and often having receivedmedications. We cannot be sure that the changes wehave identified are direct effects of the illness or are sec-ondary to an underlying pathology. We were unable toobtain medication information for all samples used inthis study (specifically, GSE17612 and Mclean66), andtherefore were unable to precisely identify the extent towhich antipsychotic use affects the results of our net-work analysis. We found that medication-associatedgenes overlap to a similar degree with the immune re-sponse module in both networks, which suggests thatcoexpression clustering patterns in schizophrenia arenot driven by medication effects. The effects of otherconfounding variables (i.e. age, pH and sex) were alsoaddressed in a similar manner, however we cannot ex-clude with certitude the possibility that the network prop-erties we have identified are still in some way influenced bythese extraneous factors. Also, samples used in our studycomprise a heterogeneous collection of cell types from theDLPFC. While the majority of our datasets utilized samplesfrom Brodmann areas 10 and 46, we also included onedataset from Brodmann area 9. Focusing on samples withina specific region would be ideal, as we avoid the potentialdilution of cell-specific biological signals associated withschizophrenia. However, we included samples from all threeBrodmann areas to maximize total sample size in our studyand increase the power of our analysis.Finally, we wish to stress limitations to the interpret-ation of coexpression networks. In contrast with otherbiological networks (i.e. protein interaction networks ormetabolic networks) whose edges represent well-definedbiological interactions, the edges in a coexpression net-work are a reduced representation of the correlationstructure of the data. The edges are related to values ofthe pairwise correlation coefficient that are calculatedfrom the expression data of the genes, and are dependenton the threshold applied to infer those networks. A connec-tion between two genes in a coexpression network does notnecessarily correspond with a connection in PPI networks,pathway or regulatory networks [46]. Thus, when studyinggene coexpression networks it is important not to confusethe edges as direct physical interactions. Indeed, there isgrowing evidence that a substantial proportion of the vari-ance in gene expression among replicate brain samples isexplained by variance in cellular composition [5,6]. Ourwork supports this as we identify “modules” that arestrongly associated with different cell type marker genes.Thus rather than identifying “physical interactions” amonggene products, coexpression patterns to a large degree seemto reflect cell-type enriched expression, that is, expressionin the same cell, but not necessarily finer levels of granular-ity such as a pathway. Changes in the composition of suchmodules might reflect changes in cellular states in subpop-ulations of cells, or changes in the associations of cell typeswith one another [47]. It is important not to interpret suchdifferences as meaning that a physical interaction has beengained or lost among gene products.ConclusionsIn summary we have contributed the largest meta-analysis of gene coexpression in schizophrenia. We eval-uated various topological properties of the control andschizophrenia networks to reveal a shared coexpressionstructure between them. Characterization of functionalclusters in each network with cell-type marker genesdisplayed differences that link together disease-relatedMistry et al. BMC Neuroscience 2013, 14:105 Page 11 of 16http://www.biomedcentral.com/1471-2202/14/105processes. Differentially expressed genes in schizophre-nia also associate with biologically relevant clustersproviding evidence for systems level dysfunction. Furtherresearch is required to disentangle these networkfindings to distinguish primary from secondary diseasephenomena, but we hope our study will encourage newdirections in the network biology of schizophrenia. Finally,our work demonstrates novel methodological approachesthat can assist in ensuring that coexpression analyses yieldbiologically justifiable and robust results.MethodsData processing and quality controlExpression profiling data sets were selected on the basisof microarray platform, use of prefrontal cortex (BA 9,10 or 46), the availability of information on covariatessuch as age, and finally the availability of the raw data[2]. Details on each of the seven datasets, including thesource citation, can be found in Table 6. Data werepreprocessed as described; briefly, expression levels weresummarized, log2 transformed and normalized for eachindividual dataset using the R Bioconductor ‘affy’ pack-age [48], with default settings for the RMA algorithm.Sample outliers were removed from each dataset basedon an inter-sample correlation analysis, resulting in atotal of 306 samples (153 from schizophrenia subjects,153 from unaffected controls) across the seven data sets.Schizophrenia and control groups had no significantdifferences in age and PMI (Table 1). Sex differenceswere assessed by use of a chi-squared test for equality ofproportion, and we observed no significant difference(p = 0.1). Brain pH was significantly different (t-test;p = 0.001). For each of the seven data sets, batch infor-mation was obtained using the ‘scan date’ stored in theCEL files; chips run on different days were considereddifferent batches and batch effects for each dataset wereremoved using the ComBat algorithm [49].Gene coexpression networksFor each dataset, samples were separated into controland schizophrenia cohorts. Probes were mapped togenes using annotations provided in Gemma (http://www.chibi.ubc.ca/Gemma), which are based on stringentmethods described in [50]. For genes mapping to mul-tiple probes, the average expression value was retained.Only genes that were represented in all seven datasetswere considered, leaving a total of 12,582 genes. Thisyielded seven expression data matrices for schizophreniaand seven for controls (one for each study). Separatenetworks were constructed for the schizophrenia andcontrol groups based on previously described methods[44]. Briefly, a gene expression profile similarity matrixwas computed for each cohort by taking the absolutevalue of the Pearson correlation between all possiblegene pairs. Correlation values in the similarity matrixwere replaced by ranks. These similarity matrices wereaggregated by cohort across datasets by taking the meanrank for each gene pair. We previously showed that thisaggregation procedure is a robust method for producinghigh-quality coexpression networks [44]. In keepingwith previous work [44,51], the aggregated matrix wasthresholded at 0.5% sparsity, resulting in an adjacencymatrix of 392,606 connections for each of the controland schizophrenia cohorts.Random coexpression networksTo evaluate the significance of network measures acrossthe whole network, formulation of appropriately ran-domized null models are required. We devised a proced-ure that results in a random network with the samenumber of genes and the same node degree distributionas the original data. Additionally, the node degree foreach individual gene is preserved (i.e. each gene still hasthe same number of connections, but the specific genesto which it is connected to are scrambled). It has beenpreviously shown that both PPI and coexpression net-works show a correlation between node degree and genemultifunctionality [45]. Thus, by constraining each geneby its node degree we can systematically assess the sig-nificance of other topological properties of the network(i.e. clustering coefficient and shortest path length) whilecontrolling for any potential multifunctionality bias inthe microarray data. To create random networks, allgene pairs were assembled into an adjacency list (2 col-umns, 392,606 rows) and genes on one side of the edgewere permuted. The resulting edges that created self-connections and/or duplicate gene pairs were isolatedand permutation was re-applied to them. This was doneiteratively until the number of conflicts was reduced toten or less. These remaining conflicting edges were re-moved from the final random network.Network propertiesWe explored three different network properties, each ofwhich is briefly described below.Node degreeEach gene can be characterized by the number of con-nections it has, that is, the number of other genes it issignificantly coexpressed with. This property is calledthe node degree. Node degrees were characterized bytheir distribution. For many biological networks the de-gree distribution has been characterized as ‘scale-free’, orat least ‘heavy tailed’. This can be observed by the qual-ity of a linear fit of the distribution on log-log scale [52].Mistry et al. BMC Neuroscience 2013, 14:105 Page 12 of 16http://www.biomedcentral.com/1471-2202/14/105Shortest pathThe shortest path length measures the shortest distanceto get from one gene to another gene by traversing edgesin the network. In an un-weighted network this is theleast number of edges traversed to get between the twogenes. We computed shortest paths using Dijkstra’s algo-rithm [53]. A value is obtained for a gene against everyother gene in the network, and presented as the meanshortest path length across all genes. Genes without anydirect neighbours are treated as missing values.Clustering coefficientThe clustering coefficient of a gene indicates howconnected the direct neighbours of a gene are to one an-other. It is the ratio of the number of connections in theneighbourhood of a node to the number of connectionsif the neighbourhood was fully connected. The clusteringcoefficient ranges from zero to one. A value of 1 would in-dicate that all the neighbours of a node are all connected toeach other, or ‘cliquish’ in nature. A value of 0 would indi-cate that none of the neighbours of a node are connectedto each other. This measure can only be computed fornodes that interact with more than one other node.Schizophrenia meta-signature network analysisThe meta-signature gene set of 25 up- and 73 down-regulated schizophrenia genes were obtained from theresults of our meta-analysis of differential expression onthe same datasets [2]. Four genes were removed fromthe down-regulated gene set as they were not present inthe network, leaving a total of 94 ‘schizophrenia genes’.Throughout this chapter we will refer to these gene setsas SZUP and SZDOWN for the genes up- and down-regulated, respectively. Average values of shortest pathlength and clustering coefficient for the SZUP andSZDOWN gene sets were evaluated within each net-work. To estimate the relevance of the network mea-sures for SZUP and SZDOWN, we implemented threeimportant controls described below.Random gene set comparisonFor each meta-signature gene set, the average values ofshortest path length and clustering coefficient werecompared to a background distribution in each network.The background distribution was generated by randomlyselecting 1000 gene sets. Random gene sets were con-strained by size and node degree of the meta-signature geneset to control for multifunctional bias. To ensure a well-matched node degree for each random gene set, selectionwas done on a per-gene basis by choosing a random genewithin ± 50 of its node degree rank. Z-scores were thencomputed to quantify the difference between the mean ofthe background distribution to the observed values for eachnetwork measure of SZUP and SZDOWN. For positivez-scores a p-value was computed reflecting how manyrandom gene sets have values higher than the observedvalue. For negative z-scores a p-value was computedreflecting how many random gene sets have values lessthan the observed value.Functional gene set comparisonAlthough our meta-signature of schizophrenia genesspan a range of cellular functions, they possess a sharedfunctional feature of altered expression in schizophrenia.Thus, it is important to assess whether the networkproperties we observe with our meta-signature gene setsare not just a property of gene groups that have lots ofshared functional features. To control for this, we gener-ated functionally characterized gene sets using the GeneOntology (GO). From the GO database (http://www.geneontology.org/), we obtained 3,230 GO terms forwhich the associated gene set size ranged from 10–1000genes. For each GO term we retrieved all human genesthat were annotated with that term to compile a genegroup for each GO term, also referred to as a functionalgene set. Each of these functional gene sets were evalu-ated individually by comparison to a background distri-bution of randomly selected gene sets (of equivalent sizeand node degree), within each network. The distributionof z-scores obtained from 3,230 functional gene groupswas plotted for each network and used to evaluate net-work properties of the meta-signature gene sets in refer-ence to other functionally related gene sets.Disease gene set comparisonTo assess the network properties of our schizophreniameta-signature genes in relation to other sets of disease-associated genes, we compiled disease gene lists for fivedifferent brain-related disorders. Gene sets were assem-bled for Alzheimer’s disease (http://www.alzgene.org/),Parkinson’s disease (http://www.pdgene.org/), multiplesclerosis (http://www.msgene.org/), and schizophrenia(http://www.szgene.org/) from their respective gene da-tabases. Each database has been compiled based on find-ings from genetic association studies and provide genelists on their website. The schizophrenia list obtainedfrom SZGene (http://www.szgene.org/) comprised onlythe top 45 of the most reliable gene associations basedon findings from a SZGene in-house meta-analysis. Wealso compiled an Autism spectrum disorder gene listfrom Toro et al. [54]. Average values of shortest pathlength and clustering coefficient were computed for allfive disease gene sets. Network measures were comparedto a background distribution of randomly selected genesets and z-scores were compared against functional geneset z-score distributions.Mistry et al. BMC Neuroscience 2013, 14:105 Page 13 of 16http://www.biomedcentral.com/1471-2202/14/105Network clusteringTo extract modules (i.e. subset of nodes that are moredensely connected to each other than to nodes outsidethe subset) from the control and schizophrenia networkswe implemented a cluster-based algorithm based onmethods described by [5]. Each adjacency matrix wastransformed into a distance matrix by computing thetopological overlap between all probe pairs [55]. Topo-logical overlap measure (TOM) between two genes iscalculated by comparing the direct connections of each.If two nodes connect to the same group of other nodesthey are said to have ‘high topological overlap’. We useda generalization of this measure that enriches TOM’ssensitivity to longer ranging connections between nodesby incorporating the number of m-step neighbours(m = 2) that a pair of node share [55]. The TOM matri-ces were subjected to WGCNA-based methods [5],whereby hierarchical clustering was applied with averagelinkage, and the resulting tree was used to define net-work modules.Enrichment analysisIn order to determine which modules were associatedwith schizophrenia we looked for enrichment of differ-entially expressed genes found from our previous meta-analysis of the same data [2]. In addition to computingoverlaps of SZUP and SZDOWN with each module, wealso looked at enrichment by utilizing the t-statistic forthe disease effect of each gene. T-statistics were enteredinto a Wilcoxon rank-sum test by module and resultantp-values were corrected for multiple testing using theBenjamini-Hochberg procedure.To examine the functional roles encoded by thesemodules we used the Gene Ontology (GO) annotationsand evaluated enrichment using the over-representationanalysis (ORA) method in ErmineJ [44,56]. The ORAmethod evaluates the genes that meet a selectioncriterion to determine if there are gene sets (GO groups)which are statistically over-represented. The ORA methodrequires the entire list of genes and their associated scoresand a score threshold must be selected. Binary scores wereused to evaluate enrichment of each module; a value of 1was assigned to genes with membership in the module anda value of 0 assigned to the rest of the genes. P-valuesfor this method are computed by using the binomialapproximation to the hypergeometric distribution, andthen corrected for multiple testing using the Benjamini-Hochberg procedure. Modules were also evaluated for CNScell type enrichment by cross-referencing the genes in eachmodule with published lists of neuron, oligodendrocyte andastrocyte marker genes [15]. Cell type marker lists werecompiled by extracting genes at a lower stringency thresh-old (as described in [8]) and a high stringency threshold(greater than 10-fold expression change). Hypergeometricprobabilities were computed to evaluate the significance ofoverlap with cell type marker lists in each module.Additional filesAdditional file 1: Table S1. Comparing brain-related disease gene setproperties to functional GO groups.Additional file 2: Figure S1. Jackknifed network measures. For eachjackknifed network (in which one dataset is removed), we computedshortest path length and clustering coefficient for SZUP and SZDOWN. Tosummarize trends observed in the jackknife analysis, we plottedclustering coefficient, shortest path length found in the CTL and SZnetworks. Results from SZUP are found in A-B, and SZDOWN in C-D. Eachline represents a different jackknifed network, with the legend indicatingwhich dataset was removed.Additional file 3: Table S2. Disease characterization of all coexpressionmodules in each network. Modules in each network were characterizedby enrichment of genes that are differentially expressed in schizophrenia.Disease effect p-values were computed by entering the t-statistic for thedisease effect of each gene into a Wilcoxon rank-sum test by module.SZUP and SZDOWN are the up- and down-regulated schizophrenia genesets previously identified in our meta-analysis [2].Additional file 4: Table S3. Cell-type marker enrichment ofcoexpression modules in each network. For each cluster in the controland schizophrenia networks we report the number of genes that overlapwith published lists of cell-type marker genes for oligodendrocytes ,neurons, and astrocyte marker genes provided by [15]. Cell-typeenrichment was computed for a low-stringency list (> 4-fold) and for ahigh stringency list (> 10-fold). Hypergeometric probabilities werecomputed to evaluate significance of overlap.Additional file 5: Table S4. Gene Ontology enrichment of top fivedisease modules in control.Additional file 6: Table S5. Gene Ontology enrichment of top fivedisease modules in schizophrenia networks, respectively.Additional file 7: Table S6. Enrichment of genes previously associatedwith other covariates.Additional file 8: Figure S2. Evaluating the effects of covariates onnetwork modules. For the age up- and pH down regulated genes whichare enriched in the immune response module of CTL, the expressiondata was plotted to evaluate differential expression between control andschizophrenia. A) Genes which remain in the SZ immune module; B)Genes that are lost from the SZ immune module. In either case, theexpression for these genes is variable within each cohort and differencesin mean expression between cohorts are very small and not significant.AbbreviationsCTL: Control network; SZ: Schizophrenia network; SZUP: Meta-signaturegenes up-regulated in schizophrenia; SZDOWN: Meta-signature genes down-regulated in schizophrenia; CNS: Central nervous system; PPI: Protein-proteininteraction; GO: Gene ontology.Competing interestsThe authors declare that they have no competing interests.Authors’ contributionsMM implemented algorithms and performed all of the analysis. JGconstructed the rank aggregated coexpression networks and designedmethods. PP oversaw design and execution of the study. PP and MM wrotethe manuscript. All authors read and approved the final manuscript.Author details1Canadian Institute of Health Research/Michael Smith Foundation for HealthResearch (CIHR/MSFHR) Graduate Program in Bioinformatics, University ofBritish Columbia, Vancouver, BC, Canada. 2Stanley Institute for CognitiveGenomics, Cold Spring Harbor Laboratory, One Bungtown Road, NY 11724,USA. 3Department of Psychiatry, University of British Columbia, Vancouver,BC, Canada. 4Centre for High-throughput Biology, 177 Michael SmithMistry et al. BMC Neuroscience 2013, 14:105 Page 14 of 16http://www.biomedcentral.com/1471-2202/14/105Laboratories, 2185 East Mall, University of British Columbia, Vancouver, BC,Canada.Received: 5 February 2013 Accepted: 23 September 2013Published: 26 September 2013References1. Sequeira PA, Martin MV, Vawter MP: The first decade and beyond oftranscriptional profiling in schizophrenia. Neurobiol Dis 2012, 45(1):23–36.2. Mistry M, Gillis J, Pavlidis P: Genome-wide expression profiling ofschizophrenia using a large combined cohort. Mol Psychiatry 2013,18(2):215–225.3. Bader GD, Hogue CW: An automated method for finding molecularcomplexes in large protein interaction networks. BMC Bioinformatics 2003, 4:2.4. Spirin V, Mirny LA: Protein complexes and functional modules inmolecular networks. Proc Natl Acad Sci U S A 2003, 100(21):12123–12128.5. Zhang B, Horvath S: A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol 2005, 4:Article17.6. Oldham MC, Konopka G, Iwamoto K, Langfelder P, Kato T, Horvath S,Geschwind DH: Functional organization of the transcriptome in humanbrain. Nat Neurosci 2008, 11(11):1271–1282.7. Gaiteri C, Sibille E: Differentially expressed genes in major depressionreside on the periphery of resilient gene coexpression networks.Frontiers in neuroscience 2011, 5:95.8. Torkamani A, Dean B, Schork NJ, Thomas EA: Coexpression networkanalysis of neural tissue reveals perturbations in developmentalprocesses in schizophrenia. Genome Res 2010, 20(4):403–412.9. Voineagu I, Wang X, Johnston P, Lowe JK, Tian Y, Horvath S, Mill J, CantorRM, Blencowe BJ, Geschwind DH: Transcriptomic analysis of autistic brainreveals convergent molecular pathology. Nature 2011, 474(7351):380–384.10. Roussos P, Katsel P, Davis KL, Siever LJ, Haroutunian V: A system-leveltranscriptomic analysis of schizophrenia using postmortem brain tissuesamples. Arch Gen Psychiatry 2012, 69(12):1–11.11. Chen C, Cheng L, Grennan K, Pibiri F, Zhang C, Badner JA, Gershon ES, Liu C,Members of the Bipolar Disorder Genome Study C: Two gene co-expressionmodules differentiate psychotics and controls. Mol Psychiatry 2012.doi:10.1038/mp.2012.146.12. Stumpf MP, Porter MA: Mathematics. Critical truths about power laws.Science 2012, 335(6069):665–666.13. Tong AH, Lesage G, Bader GD, Ding H, Xu H, Xin X, Young J, Berriz GF, BrostRL, Chang M, et al: Global mapping of the yeast genetic interactionnetwork. Science 2004, 303(5659):808–813.14. Yu H, Braun P, Yildirim MA, Lemmens I, Venkatesan K, Sahalie J, Hirozane-Kishikawa T, Gebreab F, Li N, Simonis N, et al: High-quality binary proteininteraction map of the yeast interactome network. Science 2008,322(5898):104–110.15. Cahoy JD, Emery B, Kaushal A, Foo LC, Zamanian JL, Christopherson KS,Xing Y, Lubischer JL, Krieg PA, Krupenko SA, et al: A transcriptomedatabase for astrocytes, neurons, and oligodendrocytes: a new resourcefor understanding brain development and function. J Neurosci 2008,28(1):264–278.16. Mistry M, Pavlidis P: A cross-laboratory comparison of expression profiling datafrom normal human postmortem brain. Neuroscience 2010, 167(2):384–395.17. Iwamoto K, Bundo M, Kato T: Altered expression of mitochondria-relatedgenes in postmortem brains of patients with bipolar disorder orschizophrenia, as revealed by large-scale DNA microarray analysis. HumMol Genet 2005, 14(2):241–253.18. Prabakaran S, Swatton JE, Ryan MM, Huffaker SJ, Huang JT, Griffin JL,Wayland M, Freeman T, Dudbridge F, Lilley KS, et al: Mitochondrialdysfunction in schizophrenia: evidence for compromised brainmetabolism and oxidative stress. Mol Psychiatry 2004, 9(7):684–697, 643.19. Garbett K, Gal-Chis R, Gaszner G, Lewis DA, Mirnics K: Transcriptomealterations in the prefrontal cortex of subjects with schizophrenia whocommitted suicide. Neuropsychopharmacol Hung 2008, 10(1):9–14.20. Katsel P, Davis KL, Gorman JM, Haroutunian V: Variations in differentialgene expression patterns across multiple brain regions in schizophrenia.Schizophr Res 2005, 77(2–3):241–252.21. Maycox PR, Kelly F, Taylor A, Bates S, Reid J, Logendra R, Barnes MR,Larminie C, Jones N, Lennon M, et al: Analysis of gene expression in twolarge schizophrenia cohorts identifies multiple changes associated withnerve terminal function. Mol Psychiatry 2009, 14(12):1083–1094.22. Narayan S, Tang B, Head SR, Gilmartin TJ, Sutcliffe JG, Dean B, Thomas EA:Molecular profiles of schizophrenia in the CNS at different stages ofillness. Brain research 2008, 1239:235–248.23. Walterfang M, Wood SJ, Velakoulis D, Pantelis C: Neuropathological,neurogenetic and neuroimaging evidence for white matter pathology inschizophrenia. Neurosci Biobehav Rev 2006, 30(7):918–948.24. Hansen T, Hemmingsen RP, Wang AG, Olsen L, Timm S, Soeby K, JakobsenKD, Fenger M, Parnas J, Rasmussen HB, et al: Apolipoprotein D isassociated with long-term outcome in patients with schizophrenia.Pharmacogenomics J 2006, 6(2):120–125.25. Qin W, Gao J, Xing Q, Yang J, Qian X, Li X, Guo Z, Chen H, Wang L, Huang X, etal: A family-based association study of PLP1 and schizophrenia.Neurosci Lett 2005, 375(3):207–210.26. Yang YF, Qin W, Shugart YY, He G, Liu XM, Zhou J, Zhao XZ, Chen Q, La YJ,Xu YF, et al: Possible association of the MAG locus with schizophrenia ina Chinese Han cohort of family trios. Schizophr Res 2005, 75(1):11–19.27. Bedard A, Tremblay P, Chernomoretz A, Vallieres L: Identification of genespreferentially expressed by microglia and upregulated during cuprizone-induced inflammation. Glia 2007, 55(8):777–789.28. Farina C, Aloisi F, Meinl E: Astrocytes are active players in cerebral innateimmunity. Trends Immunol 2007, 28(3):138–145.29. Condamine T, Le Texier L, Howie D, Lavault A, Hill M, Halary F, Cobbold S,Waldmann H, Cuturi MC, Chiffoleau E: Tmem176B and Tmem176A areassociated with the immature state of dendritic cells. J Leukoc Biol 2010,88(3):507–515.30. D’Agostino PM, Gottfried-Blackmore A, Anandasabapathy N, Bulloch K: Braindendritic cells: biology and pathology. Acta Neuropathol 2012, 124(5):599–614.31. Perez-Santiago J, Diez-Alarcia R, Callado LF, Zhang JX, Chana G, White CH,Glatt SJ, Tsuang MT, Everall IP, Meana JJ, et al: A combined analysis ofmicroarray gene expression studies of the human prefrontal cortexidentifies genes implicated in schizophrenia. J Psychiatr Res 2012,46(11):1464–1474.32. Li T, Li Z, Chen P, Zhao Q, Wang T, Huang K, Li J, Li Y, Liu J, Zeng Z, et al:Common variants in major histocompatibility complex region and TCF4gene are significantly associated with schizophrenia in Han Chinese.Biol Psychiatry 2010, 68(7):671–673.33. Stefansson H, Ophoff RA, Steinberg S, Andreassen OA, Cichon S, Rujescu D,Werge T, Pietilainen OP, Mors O, Mortensen PB, et al: Common variantsconferring risk of schizophrenia. Nature 2009, 460(7256):744–747.34. Lencz T, Morgan TV, Athanasiou M, Dain B, Reed CR, Kane JM, Kucherlapati R,Malhotra AK: Converging evidence for a pseudoautosomalcytokine receptor gene locus in schizophrenia. Mol Psychiatry2007, 12(6):572–580.35. Ozbey U, Tug E, Namli M: Interleukin-10 gene promoter polymorphism inpatients with schizophrenia in a region of East Turkey. World J BiolPsychiatry 2009, 10(4 Pt 2):461–468.36. Paul-Samojedny M, Owczarek A, Suchanek R, Kowalczyk M, Fila-Danilow A,Borkowska P, Kucia K, Kowalski J: Association study of interferon gamma(IFN-gamma) +874T/A gene polymorphism in patients with paranoidschizophrenia. J Mol Neurosci 2011, 43(3):309–315.37. Arion D, Unger T, Lewis DA, Levitt P, Mirnics K: Molecular evidence forincreased expression of genes related to immune and chaperonefunction in the prefrontal cortex in schizophrenia. Biol Psychiatry 2007,62(7):711–721.38. Saetre P, Emilsson L, Axelsson E, Kreuger J, Lindholm E, Jazin E:Inflammation-related genes up-regulated in schizophrenia brains.BMC Psychiatry 2007, 6(7):46.39. Shao L, Vawter MP: Shared gene expression alterations in schizophreniaand bipolar disorder. Biol Psychiatry 2008, 64(2):89–97.40. de Jong S, Boks MP, Fuller TF, Strengman E, Janson E, de Kovel CG, Ori AP,Vi N, Mulder F, Blom JD, et al: A gene co-expression network in wholeblood of schizophrenia patients is independent of antipsychotic-use andenriched for brain-expressed genes. PloS one 2012, 7(6):e39498.41. Gardiner EJ, Cairns MJ, Liu B, Beveridge NJ, Carr V, Kelly B, Scott RJ, TooneyPA: Gene expression analysis reveals schizophrenia-associateddysregulation of immune pathways in peripheral blood mononuclearcells. J Psychiatr Res 2012. doi:10.1016/j.jpsychires.2012.11.007.42. Sainz J, Mata I, Barrera J, Perez-Iglesias R, Varela I, Arranz MJ, Rodriguez MC,Crespo-Facorro B: Inflammatory and immune response genes havesignificantly altered expression in schizophrenia. Mol Psychiatry 2012.doi:10.1038/mp.2012.165.Mistry et al. BMC Neuroscience 2013, 14:105 Page 15 of 16http://www.biomedcentral.com/1471-2202/14/10543. Allen NC, Bagade S, McQueen MB, Ioannidis JP, Kavvoura FK, Khoury MJ,Tanzi RE, Bertram L: Systematic meta-analyses and field synopsis ofgenetic association studies in schizophrenia: the SzGene database.Nature genetics 2008, 40(7):827–834.44. Gillis J, Pavlidis P: The role of indirect connections in gene networks inpredicting function. Bioinformatics 2011, 27(13):1860–1866.45. Gillis J, Pavlidis P: The impact of multifunctional genes on “guilt byassociation” analysis. PloS one 2011, 6(2):e17258.46. Xulvi-Brunet R, Li H: Co-expression networks: graph properties andtopological comparisons. Bioinformatics 2010, 26(2):205–214.47. Tan PP, French L, Pavlidis P: Neuron-Enriched Gene Expression Patternsare Regionally Anti-Correlated with Oligodendrocyte-Enriched Patternsin the Adult Mouse and Human Brain. Frontiers in neuroscience 2013, 7:5.48. Team RDC: R: A Language and Environment for Statistical Computing. Vienna,Austria: In. Edited by Computing RFfS; 2011.49. Johnson WE, Li C, Rabinovic A: Adjusting batch effects in microarrayexpression data using empirical Bayes methods. Biostatistics 2007, 8(1):118–127.50. Barnes M, Freudenberg J, Thompson S, Aronow B, Pavlidis P: Experimentalcomparison and cross-validation of the Affymetrix and Illumina geneexpression analysis platforms. Nucleic Acids Res 2005, 33(18):5914–5923.51. Lee HK, Hsu AK, Sajdak J, Qin J, Pavlidis P: Coexpression analysis of humangenes across many microarray data sets. Genome Res 2004, 14(6):1085–1094.52. Albert RJ H, Barabasi AL: Internet: Diameter of the World-Wide Web.Nature 1999, 401:2.53. Dijkstra EW: A note on two problems in connexion with graphs.Numerische Mathematik 1959, 1:269–271.54. Toro R, Konyukh M, Delorme R, Leblond C, Chaste P, Fauchereau F,Coleman M, Leboyer M, Gillberg C, Bourgeron T: Key role for gene dosageand synaptic homeostasis in autism spectrum disorders. Trends Genet2010, 26(8):363–372.55. Yip AM, Horvath S: Gene network interconnectedness and thegeneralized topological overlap measure. BMC bioinformatics 2007, 8:22.56. Lee HK, Braynen W, Keshav K, Pavlidis P, Ermine J: tool for functionalanalysis of gene expression data sets. BMC Bioinformatics 2005, 6:269.doi:10.1186/1471-2202-14-105Cite this article as: Mistry et al.: Meta-analysis of gene coexpressionnetworks in the post-mortem prefrontal cortex of patients withschizophrenia and unaffected controls. BMC Neuroscience 2013 14:105.Submit your next manuscript to BioMed Centraland take full advantage of: • Convenient online submission• Thorough peer review• No space constraints or color figure charges• Immediate publication on acceptance• Inclusion in PubMed, CAS, Scopus and Google Scholar• Research which is freely available for redistributionSubmit your manuscript at www.biomedcentral.com/submitMistry et al. BMC Neuroscience 2013, 14:105 Page 16 of 16http://www.biomedcentral.com/1471-2202/14/105


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items