UBC Faculty Research and Publications

SABRE: a method for assessing the stability of gene modules in complex tissues and subject populations Shannon, Casey P; Chen, Virginia; Takhar, Mandeep; Hollander, Zsuzsanna; Balshaw, Robert; McManus, Bruce M; Tebbutt, Scott J; Sin, Don D; Ng, Raymond T Nov 14, 2016

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

Item Metadata


52383-12859_2016_Article_1319.pdf [ 1.17MB ]
JSON: 52383-1.0362020.json
JSON-LD: 52383-1.0362020-ld.json
RDF/XML (Pretty): 52383-1.0362020-rdf.xml
RDF/JSON: 52383-1.0362020-rdf.json
Turtle: 52383-1.0362020-turtle.txt
N-Triples: 52383-1.0362020-rdf-ntriples.txt
Original Record: 52383-1.0362020-source.json
Full Text

Full Text

METHODOLOGY ARTICLE Open AccessSABRE: a method for assessing the stabilityof gene modules in complex tissues andsubject populationsCasey P. Shannon1,6*†, Virginia Chen1,6†, Mandeep Takhar2, Zsuzsanna Hollander1,6, Robert Balshaw1,3,Bruce M. McManus1,4,6, Scott J. Tebbutt1,5,6, Don D. Sin5,6 and Raymond T. Ng1,2,6AbstractBackground: Gene network inference (GNI) algorithms can be used to identify sets of coordinately expressedgenes, termed network modules from whole transcriptome gene expression data. The identification of suchmodules has become a popular approach to systems biology, with important applications in translational research.Although diverse computational and statistical approaches have been devised to identify such modules, theirperformance behavior is still not fully understood, particularly in complex human tissues. Given humanheterogeneity, one important question is how the outputs of these computational methods are sensitive to theinput sample set, or stability. A related question is how this sensitivity depends on the size of the sample set. Wedescribe here the SABRE (Similarity Across Bootstrap RE-sampling) procedure for assessing the stability of genenetwork modules using a re-sampling strategy, introduce a novel criterion for identifying stable modules, anddemonstrate the utility of this approach in a clinically-relevant cohort, using two different gene network modulediscovery algorithms.Results: The stability of modules increased as sample size increased and stable modules were more likely to bereplicated in larger sets of samples. Random modules derived from permutated gene expression data wereconsistently unstable, as assessed by SABRE, and provide a useful baseline value for our proposed stability criterion.Gene module sets identified by different algorithms varied with respect to their stability, as assessed by SABRE.Finally, stable modules were more readily annotated in various curated gene set databases.Conclusions: The SABRE procedure and proposed stability criterion may provide guidance when designing systemsbiology studies in complex human disease and tissues.Keywords: Systems biology, Gene modules, Reproducibility, WGCNA, BootstrapBackgroundGene network inference (GNI) from whole transcrip-tome expression data is a fundamental and challengingtask in computational systems biology, with potentiallyimportant applications in translational research. Often, amajor aim is to identify sets of coordinately expressedgenes, termed network modules, and to study theirproperties and how they change across conditions [1–6].These gene network modules may represent novel,context specific, functional biological units, andidentifying and studying such modules has become animportant tool of systems biology [7–10]. Althoughdiverse computational and statistical approaches havebeen devised to identify such modules [11–18], theirperformance behavior is still not fully understood,particularly in complex human tissues. This situationis understandable – objective assessment is challen-ging in this context. Nevertheless, these tools arebeing applied to the study of complex diseases [19, 20],and evaluating the performance, as well as understandingthe limitations of state-of-the-art methods in this context,* Correspondence: casey.shannon@hli.ubc.ca†Equal contributors1PROOF Centre of Excellence, Vancouver, BC, Canada6UBC James Hogg Research Centre and Institute of HEART + LUNG Health,Vancouver, BC, CanadaFull list of author information is available at the end of the article© The Author(s). 2016 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link tothe Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.Shannon et al. BMC Bioinformatics  (2016) 17:460 DOI 10.1186/s12859-016-1319-8is important. This is particularly true for translational re-search where studies tend to be smaller, and are often per-formed in complex tissues and/or highly heterogeneouspatient populations. While it is challenging to directly as-sess the accuracy of inferred networks derived by thesemethods in real-world expression datasets, the resultinggene modules have properties that may be objectivelymeasured in such data.One such property is module reproducibility. A com-prehensive review of various means of assessing modulereproducibility was published by Langfelder et al. [21].In that study, the authors outlined two broad categoriesof module preservation statistics: cross-tabulation andnetwork topology derived (including density, connectiv-ity, and separability), and chose to focus on the latter asa means of assessing module reproducibility, in partbecause of the difficulty of assessing cross-tabulationmodule preservation statistics in the absence of an inde-pendent validation dataset. Re-sampling approaches,such as cross-validation or bootstrapping, have previouslybeen used to overcome this, however [21–23]. Cross-tabulation reproducibility measures are of particular inter-est because they can be readily applied to a broad range ofgene module identification strategies, including popularclustering approaches [24–26] that do not yield a nodes-and-edges network structure and, therefore, cannot beassessed by network topological metrics.Here, we introduce a flexible strategy for assessing thereproducibility of gene modules that does not rely onnetwork topology and leverages bootstrap re-sampling inplace of an independent validation dataset. Even thoughthe idea of using bootstrap re-sampling to evaluate mod-ule stability is not new (e.g., [21–23]), a systematic ap-proach to summarize the results obtained across manyre-samplings into an easy-to-interpret criterion ofmodule stability, has yet to be described. We demon-strate that the proposed procedure can provide a usefulmeasure of module stability – how sensitive a module’sgene membership is to changes in the set of samples thatwere used to identify it – in a large (n > 200), highlyrelevant clinical dataset: gene expression profiling ofperipheral whole blood, a highly heterogeneous tissue, inchronic obstructive pulmonary disease (COPD), acomplex human disease. COPD is a progressive disease,characterized by non-reversible loss of lung function andsporadic worsening of symptoms (shortness of breath,cough, etc.) termed acute exacerbations of COPD(AECOPD). These exacerbations lead to substantialmorbidity and mortality [27]. Additionally, some pa-tients experience exacerbation episodes more fre-quently, and there is interest in understanding themolecular basis, if any, of this phenotype. Since mostexacerbations are associated with bacterial or viral re-spiratory tract infections [28], immune systemmonitoring using whole blood, transcriptome-wide,gene expression profiling could provide useful insights[25]. We explored this notion by identifying gene net-work modules from this gene expression data usingthree popular strategies: weighted gene co-expressionanalysis (WGCNA) [11], a method based on the par-tial least squares regression (PLS) technique describedby Pihur, Datta and Datta [18], and a clustering-basedapproach described by Chaussabel et al. [24].The proposed procedure is more formally defined below(cf. Methods section). Briefly, the ability to recover similarmodules across many re-sampled datasets reflects modulestability. We propose to estimate gene module stability bylooking for concordance between a reference module setderived from the entire data and a large number of com-parator module sets derived from re-sampled datasets.Concordance can be determined using some similaritymeasure; we propose a variation on the Jaccard similaritycoefficient, a statistic commonly used to assess the similar-ity of sets. The stability of a particular module can be esti-mated by inspecting the distribution of similarity scoresacross a large number of repeated re-samplings. A sche-matic representation of the procedure, termed SABRE(Similarity Across Bootstrap RE-samplings) is shown inAdditional file 1: Figure S1.To facilitate the prioritizing of modules, we introduce acriterion that summarizes the distribution of similaritymeasures obtained for a particular reference moduleacross bootstrap re-samplings, and explore some usefulproperties of this criterion. We hypothesized that modulesidentified by various algorithms, such as WGCNA and theapproaches described by Pihur, Datta and Datta, orChaussabel and others, vary with respect to their stabilityand that this quality may provide a useful means of priori-tizing specific modules for further study. Moreover,changes in module stability across conditions may suggestloss of regulation, which could be of biological interest.MethodsDatasetWe obtained PAXgene blood samples from 238patients with chronic obstructive pulmonary disease(COPD), who were enrolled in the Evaluation of COPDLongitudinally to Identify Predictive Surrogate Endpoints(ECLIPSE) study [29]. These patients were clinically stableat the time of blood collection and their demographics aresummarized in Additional file 2: Table S1.RNA extraction and microarray processingBlood samples for all subjects and timepoints werecollected in PAXgene tubes and stored at −80 °C untilanalysis. Total RNA was extracted using PAXgene BloodRNA Kits (QIAGEN Inc., Germantown, MD, USA), andintegrity and concentration determined using an AgilentShannon et al. BMC Bioinformatics  (2016) 17:460 Page 2 of 112100 BioAnalyzer (Agilent Technologies Inc., SantaClara, CA, USA). Affymetrix Human Gene 1.1 ST(Affymetrix, Inc., Santa Clara, CA, USA) microarrayswere processed at the Scripps Research InstituteMicroarray Core Facility (San Diego, CA, USA) in orderto assess whole transcriptome expression. The microar-rays were checked for quality using the RMAExpresssoftware [30] (v1.1.0). All microarrays that passed qualitycontrol were background corrected and normalizedusing quantile normalization (as in RMA) [30] andsummarized using a factor analysis model (factoranalysis for robust microarray summarization [FARMS])[31], via the ‘farms’ R package. FARMS includes anobjective feature filtering technique that uses themultiple probes measuring the same target transcript asrepeated measures to quantify the signal-to-noise ratioof that specific probe set. Informative probe sets, asidentified by FARMS (2512), were used for all down-stream analyses. Limiting the feature space in thismanner had the additional benefit of speeding up identi-fication of gene network modules in the next step, animportant consideration given the proposed bootstrap-ping procedure.Identification of gene modulesWe identified network modules using three different ap-proaches: weighted gene co-expression network analysis(WGCNA) [11], via the ‘WGCNA’ R package [32], amethod based on the partial least squares regression(PLS) technique described by Pihur, Datta and Datta[18], and a k-means clustering-based approach toidentifying sets of coordinately expressed transcriptsdescribed by Chaussabel and others [24]. In WGCNA,the unsigned gene-gene correlation matrix was weightedto produce a network with approximately scale-freetopology, characteristic of biological systems [32].Average linkage hierarchical clustering is then applied tothis weighted gene co-expression network, and theresulting dendrogram cut at a height of 0.15 to producemodules with minimum size of 50. The total number ofmodules identified in this manner was not fixed. Weadopted a similar strategy to identify modules using theapproach described by Pihur and others, this timeapplying average linkage clustering to the gene-geneinteraction matrix output from the PLS procedure andcutting the resulting dendrogram at a heaight of 0.15, asbefore. The Chaussabel approach employs the k-meansclustering algorithm to identify co-clustering genes. Thisalgorithm uses a top down strategy in which genes arerandomly divided into a predetermined number ofclusters. Genes are then iteratively re-assigned to theirnearest cluster by some distance function (Euclideandistance, in this case), and cluster centers re-computedunder the new configuration. This process is repeated untilthe algorithm converges. The number of clusters, or genemodules, was necessarily fixed in this case. We used theelbow criterion to determine the point at which the inclu-sion of additional clusters no longer provides a large in-crease in proportion of variance explained. Because theinitial centers are arbitrarily chosen, the solution is unlikelyto be the minimal sum of squares of all possible partitions.Rather a local minimum is returned where any further re-assignment of a gene from one cluster to another will notreduce the within cluster sum of squares. To address thiswe used 10 random initial configurations and retained thesolution with minimal sum of squares. Although this pro-vides a minimum over several partitions, it does not guar-antee a global minimum solution.Assessing the stability of gene modulesAll three approaches were applied as described to all 238peripheral whole blood expression profiles to identifymodules of coordinately expressed transcripts. These aretermed the reference module sets for each approach. Toassess stability of all modules within the referencemodule sets, bootstrap re-sampling was used to generate1000 random sets of 238 subjects. All three networkmodule identification approaches were applied to allre-samplings to produce 1000 bootstrapped sets ofmodules, in each case.In order to describe how we propose to assess the sta-bility of the reference modules, we must first introduce afew related concepts. We start with the concept ofmodule accuracy, defined as follows: for each referencemodule q, a set of genes of size |q|, and a test module q′, where |q ∩ q′| is the number of genes in commonbetween the two. Accuracy is then defined as:Accuracyq ¼q∩q′j jjqj ∈ 0; 1½  ð1ÞWe instead propose to use the closely related Jaccardsimilarity coefficient, a statistic commonly used to assessthe similarity of sets:Jaccardqq′ ¼ q∩q′j jjq∪q′j ∈ 0; 1½  ð2ÞComparison of a reference module to a test modulethat is an exact subset should be considered a perfectmatch from a stability standpoint (the similarity measureused should return 1 in this case). In order for thesimilarity measure used to reflect this, we modify (2) asfollows:Similarityqq′ ¼q∩q′j jminfjqj; jq′jg ∈ 0; 1½  ð3ÞThis is analogous to the Simpson index that can becomputed for bipartite networks [33], but note thatShannon et al. BMC Bioinformatics  (2016) 17:460 Page 3 of 11similarity is defined only in terms of the degree ofoverlap between modules members or nodes. The mod-ules are treated as simple sets in (3).Finally, gene module stability can be very naturallyestimated via a random re-sampling with replacementprocedure, termed bootstrapping [34], by looking forconcordance between a reference module set derivedfrom the entire data and one derived from a re-sampled dataset. Concordance can be determinedusing the similarity measure defined above (3). Whencomparing a module q to a module set Q = {q1, q2,q3,…,qn}, the best match for q is defined to be q′,such that (3) between q and q′ is the highest amongall possible comparisons between q and members ofQ. The best match similarity coefficient is then (3)between q and q′, and the stability of a particularmodule q can be estimated by looking at the distribu-tion of best match similarity scores across a largenumber of repeated re-samplings, {Rj; j = 1,2,…,n}. Inorder to rank modules by their stability, we summarizethese distributions for each module to a single number,the Hirsch index (H-index), as follows:H−index qð Þ ¼ maxh 11000X1000j¼1 Rj ≥ h  ≥ h ∈ 0; 1½ ð4ÞFor a reference module with H-index = 0.8, 0.8 similar-ity or greater was observed in 80 % of bootstrap runs. Amore qualitative interpretation would be that we expectthis reference module, derived from all available sam-ples, to have 80 % similarity to a hypothetical modulederived from a whole population dataset.Construction of random gene modulesTo get a sense of the stability that could be expectedof a module containing genes with minimal relationto each other, we carried out a simulation study.Modules of size 50–400 (by increments of 50) werecreated by sampling from the all 2512 gene symbolsin the FARMS filtered dataset. This was done 100times for each size of module. Then, for each of theserandom modules, the best match Jaccard similaritycoefficient was recorded across all 1000 module setsgenerated during the previously described bootstrapprocedure. The resulting distribution was summarizedusing the H-index.Stability of network modules, sample size, and modulesizeIn order to study the relationship between networkmodule stability, sample size, and module size, a slightvariation of the bootstrap re-sampling strategy was used.We randomly sampled, without replacement 10, 20, 40,80, 120, or 160 expression profiles from the 238 peripheralwhole blood expression profiles described above. In eachcase, a reference module set was produced, 100 bootstrapre-samplings of the selected expression profiles generated,and the stability of the reference module set acrossbootstrap re-samplings determined as before. This wasrepeated 10 times to capture the effect the originalselection had on generation of the reference module set.Stability of network modules and network topologyThe relationship between module stability and variousnetwork topology measures was also of interest. Weconstructed an undirected network using the ‘igraph’ Rpackage [35]. We defined a gene-gene edge as thatwhere the absolute correlation for that gene pair was atleast two standard deviations away from the mean cor-relation observed across all possible gene pairs. Varioustopology measures were then calculated for each of thereference modules: average number of neighbors pergene (divided by module size), number of instances inwhich a gene appears in a shortest path between twoother genes, and number of triads (divided by number ofpossible triads in the module), and compared to theirstability.Stability of network modules and functional annotationFinally, we explored the relationship between stabilityand our ability to functionally annotate gene modules.We hypothesized that stable or reproducible gene mod-ules should correspond well to known biological pro-grams more often than unstable ones. To test this, wecompared the gene membership of our reference mod-ules to the MSigDB collections (v5.0) [36]. We also in-cluded the recently described Blood TranscriptomeModules (BTMs), a collection of blood-specific tran-scriptomic modules derived from an analysis of over30,000 human blood transcriptome profiles from morethan 500 studies whose data are publicly available [26].Annotation was done by testing for over-representation ofgenes from the MSigDB gene sets using a hypergeometrictest, a simple statistical method commonly used to quanti-tatively measure enrichment [37]. To quantify how well aparticular module was annotated in the gene set collec-tion, we computed the sum of the –log10 of the p-valuesfor the hypergeometric test across all gene sets in the col-lection and compared it to its stability, separately in eachof the collections (MSigDB Hallmark, C1-C7, and BTMs).ResultsWe first identified sets of gene modules by using allavailable samples and by employing three different ap-proaches: weighted gene co-expression network analysis(WGCNA) [11], a method based on the partial leastsquares regression (PLS) technique described by Pihur,Shannon et al. BMC Bioinformatics  (2016) 17:460 Page 4 of 11Datta and Datta [18], and a clustering-based approach toidentifying sets of coordinately expressed transcripts de-scribed by Chaussabel and others [24]. These are termedreference module sets.When applied to all 238 samples, WGCNA identified 19network modules, ranging in size from 69 to 850 probe-sets (36 to 427 genes; Additional file 3: Table S2), with amean module size of 240 probe-sets (median = 157). Themethod from Pihur, Datta and Datta, identified 24 networkmodules, ranging in size from 69 to 554 probe-sets (32 to210 genes; Additional file 4: Table S3), with a mean modulesize of 207 probe-sets (median = 157). The Chaussabelapproach identified 20 network modules ranging in size from44 to 302 probe-sets (38 to 167 genes; Additional file 5:Table S4), with mean module size of 131 probe-sets(median = 120). Modules identified by WGCNA and thePihur approach were highly concordant (mean Jaccardsimilarity = 0.78; not shown), while those identified by theChaussabel approach were largely distinct (Jaccard similarity< 0.2 for the majority of module pair-wise comparisons;Additional file 6: Figure S2) with the exception of theWGCNA-derived turquoise and blue modules, which weresimilar to the Chaussabel-derived M17 and M19 modules,respectively (similarity coefficient > 0.5). Annotation of theturquoise/M17, and blue/M19 modules, was consistent(Additional file 7: Table S5).Assessing gene module stability using SABRE and theH-INDEXNext, we applied SABRE to assess the stability of thesereference modules. Refer to Methods for detail. Briefly,each algorithm was allowed to identify a set of networkmodules from each re-sampling and these re-sampledmodule sets were compared to the reference module set.For each reference module, only the highest observedsimilarity coefficient within each re-sampling (that forthe best matching re-sampled module) was recorded andthe distribution of these similarity coefficients across allbootstrap re-samplings was used to assess stability.Modules with distribution of similarity coefficientsacross all re-samplings that skewed towards one wereconsistently matched to highly concordant re-sampledmodules across many random sample configurationsand are thus deemed to be relatively insensitive to sam-ple outliers. To simplify interpretation, we summarizedthe distribution of similarity coefficients for each moduleto its Hirsch-index [38] (H-index), as defined in (4). Avisual derivation of the H-index is shown in Fig. 1 forthe modules identified by WGCNA.Gene module stability increases with sample size andmodule sizeGene modules are often identified in relatively smallstudy populations. We used a variation of the SABREstrategy to study the effect of sample size on gene modulestability. As expected, module stability, as assessed usingthe current framework, increased as the sample size was in-creased (Fig. 2a). This relationship held for both very stable(1st rank) and less stable (10th rank) modules. We alsoobserved a relationship between module stability andmodule size in smaller studies (Fig. 2b; n= 10, p= 4.3 ×10−13; n= 80, p= 0.041; Wilcoxon’s rank-sum test), but therewas no such relationship at larger sample sizes (n = 120,p= 0.55; n= 160, p= 0.88). In all cases, modules identified(by WGCNA in this case) were more stable than modulesassembled by randomly sampling from all gene symbols inthe filtered dataset, even when sample size was small(Additional file 8: Figure S3).Stability profiles differ between algorithmsThe H-index was next used to rank reference modules(Fig. 3). Modules varied with respect to their stability asassessed by the SABRE procedure, both within sets ofmodules identified by particular algorithms, as well asacross strategies. Modules identified by all three strat-egies were significantly more stable than modules assem-bled by randomly sampling from all gene symbols in thefiltered dataset (Additional file 8: Figure S3). The set ofnetwork modules identified by WGCNA was morestable (mean H-index = 0.79, standard deviation = 0.12,range = 0.42–0.96) than that identified by either thePihur (mean H-index = 0.68, standard deviation = 0.12,range = 0.44–0.91), or Chaussabel approaches (meanH-index = 0.69, standard deviation = 0.14, range = 0.43–0.97), though the top ranking module identified by allthree approaches had comparable stability (WGCNA:lightgreen module, H-index = 0.96; Pihur: blue module,H-index = 0.91; Chaussabel: M6 module, H-index = 0.97).Stable modules are more interconnectedWe might expect more connected gene network mod-ules to be more stable. There are many proposed metricsfor measuring network connectivity. The simplest suchmetric is the average number of neighbors. For our data-set, there is no observed trend between the H-index andthe average number of neighbors within a module. Astronger notion of connectivity is to measure the num-ber of nodes that are in the shortest path between somepair of nodes in the network. The higher this number,the more tightly connected the network. Figure 4 sug-gests that there is a relationship between this notion ofconnectivity and the H-index. An even stronger notionof connectivity described in the literature is the numberof triads, which are triplets of nodes that are directlyconnected pairwise. Figure 4 indicates that there is aweak relationship between this notion of connectivityand the H-index, particularly for lower values of the H-index. In sum, it appears that the notion of stabilityShannon et al. BMC Bioinformatics  (2016) 17:460 Page 5 of 11indicated by the H-index is related, but not identical, to 2well-known notions of connectivity in network science.Stable modules are more readily annotatedWe hypothesized above that stable modules should cor-respond to well-characterized biological functions. If thisis true, we would expect stable modules to be morereadily annotated than less stable ones. We explored thisnotion in the BTM and MSigDB collections of annotatedgene sets and found that stable modules were indeedmore readily annotated across many of the includedcollections (Fig. 5). Module stability was significantlyassociated with annotatability (sum –log10 p-value ofthe hypergeometric test of the overlap betweenannotation and module genes) in the hallmark (H,Spearman’s ρ = 0.39; p = 0.01), positional (C1, Spear-man’s ρ = 0.39; p = 0.01), immunologic signatures (C7,Spearman’s ρ = 0.31; p = 0.05), and blood transcriptomicmodules (BTM, Spearman’s ρ = 0.31; p = 0.05) collec-tions, marginally associated in the canonical (C2,Spearman’s ρ = 0.27; p = 0.09), and oncogenic signa-tures (C6, Spearman’s ρ = 0.27; p = 0.10) collections,and showed no significant association in the geneontology (GO, C5), motif (C3), or computational (C4)gene set collections. Complete gene set over-representation results are tabulated in Additional file 5:Table S4.Many gene modules identified by WGCNA had veryhigh concordance with gene sets corresponding to dis-tinct biological functions (e.g. lightgreen: B cell activity,antibody production; cyan: interferon signaling; brown:heme metabolism; blue: recruitment of neutrophils andTLR mediated inflammatory signaling). In fact, evenrelatively less stable modules identified by WGCNA ap-peared to correspond to known biological pathways (e.g.purple module [H-index = 0.73]: MHC class II antigenpresentation). Modules identified by the Chaussabel ap-proach had generally lower concordance to annotatedgene sets, with modules often having seemingly overlap-ping biological function (e.g. interferon signaling pathwayFig. 1 Visual derivation of the H-index. A visual depiction of the derivation of the H-index is shown for modules identified by the WGCNA algorithm.A set of reference modules derived from all available samples are compared to a series of comparator module sets derived from bootstrapped re-sampled data. For each re-sampled dataset, all reference modules are compared to all newly identified modules using the Jaccard similarity coefficient.For each reference module, the best match Jaccard similarity coefficient value is recorded. Finally, these best match similarity coefficients are sortedand a measure of the area under the resulting curve (Hirsch-index) used to estimate the reference module’s stabilityShannon et al. BMC Bioinformatics  (2016) 17:460 Page 6 of 11activity was assigned to both modules M04 and M06). Infact, nearly half of gene modules identified by theChaussabel approach were enriched in genes associ-ated with translation (7/20 modules had significantenrichment with KEGG ribosome, KEGG translation) andmRNA transcription (4/20 modules had significant en-richment with KEGG spliceosome, REACTOME metabol-ism of mRNA, REACTOME Translation, REACTOMEmRNA splicing).DiscussionThe main objective of this work was to devise a strategyfor evaluating the stability of gene modules in a mannerapplicable to a broad range of gene module identificationstrategies, and not reliant on the availability of additionaldatasets for validation. Quantifying network modulestability would allow for prioritization of modules forfurther experimental interrogation. While the use of re-sampling strategies, such as cross-validation or thebootstrap, in this context has been previously proposed[21–23], a systematic approach to summarizing the re-sults obtained across many re-samplings into an easy-to-interpret criterion of module stability, has yet to be de-scribed. As shown in Fig. 1, for a particular module, allthe re-sampling results can be captured by a curve. Thequestion we ask is how to summarize the curve with asingle value that is informative. Standard ways, such asusing the mean value, or even using the area-under-the-curve value, do not provide a “minimum quality”guarantee. In contrast, the proposed bootstrap re-samplingand summarization scheme, SABRE, does provide such alower bound guarantee, at least across all generatedFig. 2 Gene module stability increases with sample size and module size. For each n, we sampled without replacement from all available geneexpression profiles 10 times. In each case, a reference module set was produced (by WCGNA), 100 bootstrap re-samplings of the selected expres-sion profiles generated, and the stability of the reference module set across bootstrap re-samplings determined as described in the Methods section.a Stability of the modules is visualized at n = 10, 20, 40, 80, 120 and 160 for the 1st, 5th and 10th rank modules. b Stability is plotted against modulesize at n = 10, 20, 40, 80, 120 and 160. The dotted line depicts the best-case stability of random modules in simulation. We compare the stability of S(1st quartile) and XL (4th quartile) modules using Wilcoxon’s rank-sum testFig. 3 Stability profiles differ between algorithms. Gene module similarity across bootstrap re-samplings, for all reference network modules identifiedby three gene module discovery algorithms, is visualized using box plots (a: WGCNA; b: Pihur, c: Chaussabel). The stability of the network modules issummarized using the H-index (red). The dotted line depicts the best-case stability of random modules in simulationShannon et al. BMC Bioinformatics  (2016) 17:460 Page 7 of 11bootstrap runs. The proposed stability criterion, theH-index of the curve, corresponding to the largestsquare under the curve, is readily interpreted: for amodule with H-index = 0.8, one can say that similarityof 0.8 or greater was observed across at least 80 % ofthe bootstrap runs.We go on to explore a number of useful characteristicsof the H-index. We show that the H-index of a modulegenerally increases as sample size increases and observethat, for any given module, a stability maximum appearsto be reached between n = 80–120, at least in this tissueand patient population. We also note that, while smallermodules are generally less stable, this is not the casewith larger sample sizes (n > 80–120). As expected, wefind randomly assembled modules to have very low sta-bility (H-index = 0.25), and it is comforting to note that,even for very small studies (n = 10), modules identifiedby WGCNA had worst case stability that exceeded thisvalue, suggesting that, for some gene modules at least,core members may be identifiable even in very smallstudies. Taken together, these observations suggest thatthe identification of robust gene expression modules incomplex tissues and diseases requires large study popu-lations (n > 100).Next, we compare gene network modules identifiedby three popular gene network module identificationFig. 4 Stable modules are more interconnected. The relationship between module stability and a number of topological measures of networkconnectivity is visualized for modules identified by WGCNA (blue) or the Chaussabel approach (red). Stability is positively associated (Spearman’sρ) with both number of appearance in the shortest path and number of triads in the network (* p≤ 0.05)Fig. 5 Stable modules are more readily annotated. Module gene over-representation in annotated gene sets (sum of –log10 p-value for thehypergeometric test) is visualized, for modules with varying stability, in the MSigDB and BTM collections. Stability is positively associated(Spearman’s ρ) with our ability to assign module to known biology (* p ≤ 0.05; † p ≤ 0.10) in many of these collectionsShannon et al. BMC Bioinformatics  (2016) 17:460 Page 8 of 11strategies. WGCNA and the PLS-based approach de-scribed by Pihur, Datta and Datta identified largelyconcordant sets of gene network modules, while themodules identified by the Chaussabel approach werelargely distinct. Given that module identification forboth WGCNA and the Pihur approach utilized aver-age linkage hierarchical clustering on an adjacencymatrix, this is perhaps not surprising. Applying SABRE tothese modules sets allowed us to readily rank identifiedgene modules from most to least stable. Modules identi-fied by WGCNA were generally more stable than thoseidentified by the other two approaches. The ranking ofmodules by their stability, in combinations with othermetrics, such as biological annotation, may informprioritization of certain modules for further study.We found that the H-index was positively associatedwith topological notions of network connectivity, as wellas our ability to assign biological function to gene mod-ules. SABRE uncovered important qualitative differencesbetween module sets in this respect, however. First,while stability was positively correlated with connectivityacross all modules, this relationship was strongest inmodules identified by the Chaussabel approach. Thesemodules also exhibited lower network connectivity com-pared to those found by WGCNA. This is not surprisingsince the Chaussabel approach pays no attention to thetopology of the constructed modules. A similar patternemerges when comparing the relationship between mod-ule stability and annotatability in the MSigDB and BTMgene set collections. Here again, the relationship be-tween stability and annotatability was strongest in mod-ules identified by the Chaussabel approach. Given thedistribution of the stability criterion, and the connectiv-ity or annotatability measures considered, it is difficultto determine whether the observed relationships be-tween stability and connectivity/annotatability are trueor primarily driven by broad differences in stability be-tween module sets identified by different strategies.Stable gene modules that do not readily correspond toannotated gene sets may be very interesting, of course,as they may represent novel, disease- or tissue-specificbiological processes. Two such gene modules were iden-tified in the current study: the lightgreen (h = 0.96) andblue (h = 0.86) modules. In an effort to assign biologicalfunction to these modules, we compared them to the re-cently published Blood Transcriptome Modules (BTM).These empirically derived sets of co-regulated geneshave very low overlap with presently available pathwaysand were identified in peripheral whole blood geneexpression data. We reasoned that the highly stable,un-annotated modules we independently identified inthis study may in fact correspond to some of thesehighly stable blood modules. In fact, this was the case:the lightgreen module was enriched for a number ofBTMs related to B cell activity, while the blue modulewas matched to various innate immunity BTMs (recruit-ment of neutrophils and TLR mediated inflammatorysignaling). Both B cells and neutrophils are known to beimplicated in COPD. [39, 40] This provided validation,both of the modules themselves, and of the stabilityranking produced by the SABRE procedure, in that twohighly stable modules that did not correspond to anyavailable pathway annotations, were consistent with in-dependently derived functional modules specific toblood leukocyte sub-populations. These modules mayrepresent important and novel biological function in theperipheral whole blood compartment of COPD patients.ConclusionsIn conclusion, we demonstrate that bootstrap re-sampling, and the SABRE procedure described herein,can assess the stability of gene modules identified bythree different algorithms and suggest that this could bea useful criterion when selecting modules for furtherinvestigation. We also show that when modules areidentified in smaller studies, more stable ones are morelikely to replicate in larger experiments compared to lessstable ones. We show a relationship between this notionof stability, topological connectivity, and our ability toassign biological function to gene modules. Our ap-proach highlights the relative robustness of the WGCNAalgorithm to sample outliers and, more generally, sug-gests that many gene module strategies should probablybe applied jointly to any given dataset. Finally, we iden-tify and validate two highly stable modules that may rep-resent novel, tissue-specific biological function in thecontext of the peripheral whole blood of clinically stableCOPD patients.Additional filesAdditional file 1: Figure S1. Schematic of bootstrap re-samplingprocedure. (PNG 146 kb)Additional file 2: Table S1. Patient demographics. Demographiccharacteristics of the 238 COPD patients. (CSV 495 bytes)Additional file 3: Table S2. WGCNA reference module sets.(CSV 550 bytes)Additional file 4: Table S3. Pihur reference module sets. (CSV 701 bytes)Additional file 5: Table S4. Chaussabel reference module sets.(CSV 497 bytes)Additional file 6: Figure S2. WGCNA and Chaussabel moduleconcordance. Concordance between network modules identified by twogene network module discovery. Reference module sets were identifiedusing all available gene expression profiles. The similarity coefficient foreach pair-wise comparison is visualized as a clustered heatmap, with redindicating high similarity and blue indicating low similarity. (PNG 77 kb)Additional file 7: Table S5. Pathway over-representation analysis resultsfor the WGCNA and Chaussabel modules sets. For each moduleidentification strategy and module in turn, we carried out pathwayover-representation analysis against the Broad Institute’s MSigDBShannon et al. BMC Bioinformatics  (2016) 17:460 Page 9 of 11collections by comparing the gene set and module membershipusing a simple hypergeometric test. We report both p-value and falsediscovery rate adjusting for multiple comparisons using the Benjamini-Hochberg procedure. (CSV 196 kb)Additional file 8: Figure S3. Random module stability. To get a senseof the stability that could be expected of a module containing geneswith minimal relation to each other, a simulation study was carried out.Modules of size 50, 100, 150, 200, 250, 300, 350, and 400 were randomlyassembled by sampling from the all 2512 gene symbols in the filtereddataset. This was done 100 times for each size of module. For each randommodule, their best match Jaccard similarity ceofficients were computed foreach of the 1000 bootstrap results previously generated, and the resultingdistribution was summarized using the h-index. (PNG 51 kb)AbbreviationsAECOPD: Acute exacerbation of chronic obstructive pulmonary disease;COPD: Chronic obstructive pulmonary disease; GNI: Gene network inference;SABRE: Stability across bootstrap re-samplings; WGCNA: Weighted gene co-expression analysisAcknowledgementsThe authors would like to thank the research participants whose tissuedonations made this work possible.FundingFunding was provided by Genome Canada, Genome British Columbia,Genome Quebec, the Canadian Institutes of Health Research, PROOF Centre,St. Paul’s Hospital Foundation, Providence Health Care, and the COPD ClinicalResearch Network.Availability of data and materialsThe gene expression data used in this study is available on the GeneExpression Omnibus: GSE71220.Authors’ contributionsCPS, VC, MT, ZH, and RTN conceived of, and designed this study. VC and ZHcarried out quality control and normalization of the microarray data. CPS, VC,and MT performed the statistical analyses. CPS drafted the manuscript. RBand SJT participated in the design of the study and helped draft themanuscript. BMM, DDS, and RTN conceived of, and designed the COPDECLIPSE gene expression profiling study whose samples were used for thisanalysis. All authors read and approved the final manuscript.Competing interestsThe authors declare that they have no competing interests.Consent for publicationNot applicable.Ethics approval and consent to participateThe ECLIPSE study was conducted in accordance with the Declaration ofHelsinki and good clinical practice guidelines, and has been approvedby the relevant ethics and review boards at the participating centres.Study participants gave written informed consent. Consent was given fortheir data to be shared as part of future research for the purpose offurther understanding the disease, such as this study. Patient informationwas anonymized and de-identified prior to the analysis. ECLIPSE study wasfunded by GlaxoSmithKline, under ClinicalTrials.gov identifier NCT00292552 andGSK study code SCO104960. This study, analyzing a subset of ECLIPSE subjects’mRNA samples, was approved by the University of British Columbia (UBC),under Research Ethics Board (REB) number H11-00786.Author details1PROOF Centre of Excellence, Vancouver, BC, Canada. 2Department ofComputer Science, UBC, Vancouver, BC, Canada. 3BC Centre for DiseaseControl, Vancouver, BC, Canada. 4Department of Pathology and LaboratoryMedicine, UBC, Vancouver, BC, Canada. 5Department of Medicine, Division ofRespiratory Medicine, UBC, Vancouver, BC, Canada. 6UBC James HoggResearch Centre and Institute of HEART + LUNG Health, Vancouver, BC,Canada.Received: 25 June 2016 Accepted: 3 November 2016References1. Kurata H, El-Samad H, Iwasaki R, Ohtake H, Doyle JC, Grigorova I, et al.Module-based analysis of robustness tradeoffs in the heat shock responsesystem. PLoS Comput Biol. 2006;2:e59.2. Xia K, Xue H, Dong D, Zhu S, Wang J, Zhang Q, et al. Identification of theproliferation/differentiation switch in the cellular network of multicellularorganisms. PLoS Comput Biol. 2006;2:e145.3. Ghazalpour A, Doss S, Zhang B, Wang S, Plaisier C, Castellanos R, et al.Integrating genetic and network analysis to characterize genes related tomouse weight. PLoS Genet. 2006;2:e130.4. Wang J, Zhang S, Wang Y, Chen L, Zhang X-S. Disease-aging networkreveals significant roles of aging genes in connecting genetic diseases. PLoSComput Biol. 2009;5:e1000521.5. Stone EA, Ayroles JF. Modulated modularity clustering as an exploratorytool for functional genomic inference. PLoS Genet. 2009;5:e1000479.6. Plaisier CL, Horvath S, Huertas-Vazquez A, Cruz-Bautista I, Herrera MF,Tusie-Luna T, et al. A systems genetics approach implicates USF1,FADS3, and other causal candidate genes for familial combinedhyperlipidemia. PLoS Genet. 2009;5:e1000642.7. Segal E, Shapira M, Regev A, Pe’er D, Botstein D, Koller D, et al. Modulenetworks: identifying regulatory modules and their condition-specificregulators from gene expression data. Nat Genet. 2003;34:166–76.8. Friedman N. Inferring cellular networks using probabilistic graphical models.Science. 2004;303:799–805.9. Basso K, Margolin AA, Stolovitzky G, Klein U, Dalla-Favera R, Califano A.Reverse engineering of regulatory networks in human B cells. Nat Genet.2005;37:382–90.10. Lee I, Lehner B, Crombie C, Wong W, Fraser AG, Marcotte EM. A single genenetwork accurately predicts phenotypic effects of gene perturbation inCaenorhabditis elegans. Nat Genet. 2008;40:181–8.11. Zhang B, Horvath S. A general framework for weighted gene co-expressionnetwork analysis. Stat Appl Genet Mol Biol. 2005;4:Article17.12. Huynh-Thu VA, Irrthum A, Wehenkel L, Geurts P. Inferring regulatorynetworks from expression data using tree-based methods. PLoS ONE.2010;5:e12776.13. Rajapakse JC, Mundra PA. Stability of building gene regulatory networkswith sparse autoregressive models. BMC Bioinformatics. 2011;12:S17.14. Haury A-C, Mordelet F, Vera-Licona P, Vert J-P. TIGRESS: trustful inference ofgene REgulation using stability selection. BMC Syst Biol. 2012;6:145.15. Kuffner R, Petri T, Tavakkolkhah P, Windhager L, Zimmer R. Inferring generegulatory networks by ANOVA. Bioinformatics. 2012;28:1376–82.16. Ruyssinck J, Huynh-Thu VA, Geurts P, Dhaene T, Demeester P, Saeys Y.NIMEFI: gene regulatory network inference using multiple ensemble featureimportance algorithms. PLoS ONE. 2014;9:e92709.17. Gill R, Datta S, Datta S. A statistical framework for differential networkanalysis from microarray data. BMC Bioinformatics. 2010;11:95.18. Pihur V, Datta S, Datta S. Reconstruction of genetic association networksfrom microarray data: a partial least squares approach. Bioinformatics.2008;24:561–8.19. de Jong S, Boks MPM, Fuller TF, Strengman E, Janson E, de Kovel CGF, et al.A Gene Co-Expression Network in Whole Blood of Schizophrenia Patients IsIndependent of Antipsychotic-Use and Enriched for Brain-Expressed Genes.Mazza M, editor. PLoS ONE. 2012;7:e39498.20. Van Eijk KR, de Jong S, Boks MP, Langeveld T, Colas F, Veldink JH, et al.Genetic analysis of DNA methylation and gene expression levels in wholeblood of healthy human subjects. BMC Genomics. 2012;13:636.21. Langfelder P, Luo R, Oldham MC, Horvath S. Is my network modulepreserved and reproducible? PLoS Comput Biol. 2011;7:e1001057.22. Brock G, Pihur V, Datta S, Datta S, others. clValid, an R package for clustervalidation. J. Stat. Softw. Brock Al March 2008 [Internet]. 2011 [cited 2016Sep 8]; Available from: http://cran.us.r-project.org/web/packages/clValid/vignettes/clValid.pdf23. Datta S, Datta S. Comparisons and validation of statistical clustering techniquesfor microarray gene expression data. Bioinformatics. 2003;19:459–66.24. Chaussabel D, Quinn C, Shen J, Patel P, Glaser C, Baldwin N, et al. A modularanalysis framework for blood genomics studies: application to systemiclupus erythematosus. Immunity. 2008;29:150–64.Shannon et al. BMC Bioinformatics  (2016) 17:460 Page 10 of 1125. Chaussabel D, Pascual V, Banchereau J. Assessing the human immunesystem through blood transcriptomics. BMC Biol. 2010;8:84.26. Li S, Rouphael N, Duraisingham S, Romero-Steiner S, Presnell S, Davis C,et al. Molecular signatures of antibody responses derived from a systemsbiology study of five human vaccines. Nat Immunol. 2013;15:195–204.27. Mannino DM, Buist AS. Global burden of COPD: risk factors, prevalence, andfuture trends. Lancet. 2007;370:765–73.28. Wedzicha JA, Seemungal TA. COPD exacerbations: defining their cause andprevention. Lancet. 2007;370:786–96.29. Vestbo J, Anderson W, Coxson HO, Crim C, Dawber F, Edwards L, et al.Evaluation of COPD longitudinally to identify predictive surrogateend-points (eclipse). Eur Respir J. 2008;31:869–73.30. Bolstad BM, Irizarry R, Åstrand M, Speed TP. A comparison of normalizationmethods for high density oligonucleotide array data based on variance andbias. Bioinformatics. 2003;19:185–93.31. Hochreiter S, Clevert D-A, Obermayer K. A new summarization method foraffymetrix probe level data. Bioinformatics. 2006;22:943–9.32. Langfelder P, Horvath S. WGCNA: an R package for weighted correlationnetwork analysis. BMC Bioinformatics. 2008;9:559.33. Fuxman Bass JI, Diallo A, Nelson J, Soto JM, Myers CL, Walhout AJM.Using networks to measure similarity between genes: association indexselection. Nat Methods. 2013;10:1169–76.34. Efron B. Bootstrap methods: another look at the Jackknife. Ann Stat.1979;7:1–26.35. Csardi G, Nepusz T. The igraph software package for complex networkresearch. Int J Complex Syst. 2006;1695:1–9.36. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA,et al. Gene set enrichment analysis: a knowledge-based approach forinterpreting genome-wide expression profiles. Proc Natl Acad Sci U S A.2005;102:15545–50.37. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools:paths toward the comprehensive functional analysis of large gene lists.Nucleic Acids Res. 2009;37:1–13.38. Hirsch JE. An index to quantify an individual’s scientific research output.Proc Natl Acad Sci U S A. 2005;102:16569–72.39. Campbell JD, McDonough JE, Zeskind JE, Hackett TL, Pechkovsky DV,Brandsma C-A, et al. A gene expression signature of emphysema-relatedlung destruction and its reversal by the tripeptide GHK. Genome Med.2012;4:67.40. Hoenderdos K, Condliffe A. The neutrophil in chronic obstructive pulmonarydisease. Too little, too late or too much, too soon? Am J Respir Cell MolBiol. 2013;48:531–9.•  We accept pre-submission inquiries •  Our selector tool helps you to find the most relevant journal•  We provide round the clock customer support •  Convenient online submission•  Thorough peer review•  Inclusion in PubMed and all major indexing services •  Maximum visibility for your researchSubmit your manuscript atwww.biomedcentral.com/submitSubmit your next manuscript to BioMed Central and we will help you at every step:Shannon et al. BMC Bioinformatics  (2016) 17:460 Page 11 of 11


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