Meta-analysis on blood transcriptomic studies identifies consistently coexpressed protein–protein interaction modules as robust markers of human aging

The bodily decline that occurs with advancing age strongly impacts on the prospects for future health and life expectancy. Despite the profound role of age in disease etiology, knowledge about the molecular mechanisms driving the process of aging in humans is limited. Here, we used an integrative network-based approach for combining multiple large-scale expression studies in blood (2539 individuals) with protein–protein Interaction (PPI) data for the detection of consistently coexpressed PPI modules that may reflect key processes that change throughout the course of normative aging. Module detection followed by a meta-analysis on chronological age identified fifteen consistently coexpressed PPI modules associated with chronological age, including a highly significant module (P = 3.5 × 10−38) enriched for ‘T-cell activation’ marking age-associated shifts in lymphocyte blood cell counts (R2 = 0.603; P = 1.9 × 10−10). Adjusting the analysis in the compendium for the ‘T-cell activation’ module showed five consistently coexpressed PPI modules that robustly associated with chronological age and included modules enriched for ‘Translational elongation’, ‘Cytolysis’ and ‘DNA metabolic process’. In an independent study of 3535 individuals, four of five modules consistently associated with chronological age, underpinning the robustness of the approach. We found three of five modules to be significantly enriched with aging-related genes, as defined by the GenAge database, and association with prospective survival at high ages for one of the modules including ASF1A. The hereby-detected age-associated and consistently coexpressed PPI modules therefore may provide a molecular basis for future research into mechanisms underlying human aging.


Introduction
A steadily growing life expectancy of the general western population throughout the past two centuries (Oeppen & Vaupel, 2002) has imposed the urgency for understanding the adverse effects of aging for public health and its relation to the observed large variation in healthy lifespan (Hitt et al., 1999). Age-dependent detrimental processes strongly attenuate prospects for future health, with chronological age being the major risk factor for mortality and virtually all common diseases in the western world (Wilson et al., 1998). Aging is a systemic ailment marked by a gradual metabolic decline eventually leading to a state of senescence on both the cellular and organismal level that seems to be caused by the accumulation of damage over time (Kirkwood, 1977). Despite their profound role for disease etiology, the existing knowledge concerning the molecular mechanisms driving biological aging processes in humans is limited.
Construction of consistent age-associated signatures has proven to be challenging as a multitude of gene expression studies have identified age-associated genes so far, though with limited mutual overlap (Passtoors et al., 2008;de Magalhaes et al., 2009). This inconsistency is most likely due to variable technical circumstances, small study sizes, and low signal-to-noise ratios, typically observed when analyzing the aging transcriptome. More similarity was observed at the pathway level, across tissues and even species (Partridge & Gems, 2002;Zahn et al., 2006), suggesting that the analysis of the aging transcriptome by functionally grouped gene sets is a promising alternative for the classical individual-gene analyses.
Rather than employing literature-based sets of genes sharing similar biological functions, so-called network approaches are increasingly used, which infer functional clusters of genes from the expression data itself by exploiting gene coexpression patterns hidden within the data (Zhang & Horvath, 2005). Alternatively, changes in these gene coexpression patterns that occur with age might be used for inferring a functional grouping from the data (Southworth et al., 2009). However, coexpression patterns may contain spurious gene-gene correlations (Stuart et al., 2003), which makes the use of multiple data sources simultaneously or the integration with other additional information sources on functional relationships between genes desirable.
Established modulators of aging processes in model organisms were reported to spatially cluster within networks constructed of proteinprotein interaction (PPI) data (de Magalhaes & Toussaint, 2004;Bell et al., 2009). Hence, PPI networks can be exploited for prioritizing new aging-associated genes (Witten & Bonchev, 2007;Tacutu et al., 2012) or for refining modules of coexpressed genes that are correlated during the course of aging (Xue et al., 2007). We previously demonstrated that the inference of these so-called coexpressed PPI modules has a high reproducibility across multiple expression datasets in breast cancer (van den Akker et al., 2011), and here we extend this algorithm to combine multiple gene expression datasets on aging.
Though many algorithms for network inference exist (Marbach et al., 2012), relatively little attention has gone to the problem of network inference and subsequent associations with a phenotype using multiple heterogeneous expression data sources simultaneously. Merging the expression data into a single set and using this for network inference clearly surpasses the differences in correlation structures present within each dataset. Irrespective of the type of network inference chosen, we propose to handle such heterogeneity by integrating the gene-gene similarity measures obtained across expression datasets using a suitable meta-analysis setting. Thus, in the approach described in this paper, we employ a meta-analysis for inferring a consistent gene-gene network that serves as a basis for identifying consistently coexpressed PPI modules, which are subsequently analyzed with respect to chronological age across datasets using again a meta-analysis.
To robustly characterize the changes of the blood transcriptome associated with chronological age, we have build a compendium using three large-scale transcriptomic studies (Goring et al., 2007;Emilsson et al., 2008;Inouye et al., 2010) generated in blood comprising 2539 individuals on which we applied our integrative network approach. For comparison, two types of individual-gene meta-analyses were performed as well, which in combination with an enrichment analysis yielded only broad terms for age-associated cellular processes. Application of our integrative network-based approach, yielded five consistently coexpressed PPI modules showing robust age associations and functional enrichments for 'Translational elongation', 'Cytolysis' and 'DNA metabolic process', which seem to reflect downstream mTOR signaling events or cell-cycle checkpoints. Finally, we show that four of five modules replicate in an independent cohort, and that they are enriched for known longevity-and aging-related genes and that the expression of one module associates with prospective survival at old age.

Results
The largest transcriptome compendium for normative aging To robustly characterize the changes of the blood transcriptome throughout the course of normative aging in the range of 15-94 years, we built a gene expression compendium using three large-scale transcriptomic studies performed in blood: the San Antonio Family Heart Study (SAFHS) (Goring et al., 2007), the Icelandic Family Blood (IFB) cohort (Emilsson et al., 2008)  Limited overlap of age-associated genes between studies within the compendium The most straightforward method for an integrative analysis across datasets is to first compute the age-association genes per dataset and subsequently inspect the overlap of significant results. A linear model adjusted for gender yielded between 111 (1.2%) and 1103 (12.2%) significantly age-associated genes per dataset (Bonferroni correction, a ≤ 0.05), of which 26 genes were significantly associated with age in all four datasets [ Fig. 1 and Table S1, Supporting information]. These results confirmed the high discrepancy between lists of age-associated genes previously reported in literature, even though now observed in equal or similar tissues (Passtoors et al., 2008;de Magalhaes et al., 2009). Rank-based integration of age-associated genes improves consistency between studies As repeatedly applied cutoffs across multiple heterogeneous datasets may lead to high false exclusion rates of age-associated genes, we investigated whether age-association rankings were consistently high across datasets by applying a rank integration approach (Breitling et al., 2004;de Magalhaes et al., 2009 Tables 2  and 3, respectively, and include many of the age-associated genes previously identified, like LRNN3, LEF1, and SYT11 (Hong et al., 2008;Harries et al., 2011;Passtoors et al., 2012). Results for all 9047 genes in the compendium are provided in Table S2 (Supporting information).

A novel integrative network approach for detecting consistent coexpressed PPI modules
To improve robustness against noise and increase power, we used a novel integrative network-based approach to explore functional age-associated groupings of genes. The proposed approach detects consistently coexpressed PPI modules across multiple datasets (for details see Experimental Procedures and Data S1, Supporting information). Using the four transcriptomic datasets mapped onto the PPI network, we detected a total of 162 consistently coexpressed PPI modules ranging in size from 2 to 37 genes [see Fig. S1, Supporting information for a complete overview]. The following steps in our analysis were limited to the subset of 27 coexpressed PPI modules counting at least five genes. Application of DAVID yielded significant functional enrichments for 19 of the 27 identified coexpressed PPI modules [ Table S5, Supporting information], suggesting that the applied approach grouped genes according to plausible biological functions.

Age-associated coexpressed PPI modules point toward T-cell activation
To test whether transcriptional changes of the 27 identified modules associate with chronological age, an expression profile for each module was constructed by determining the mean expression of the genes within a detected coexpressed PPI module per individual. As with the individual-gene analysis, we proceeded by computing the associations of the module expressions with age while adjusting for gender for each dataset separately. Only one module [ Fig. 2A], enriched for 'T-cell activation', was significantly associated with age in each of the four datasets of the compendium. This module A contains genes commonly employed as markers for assessing the differentiation status of T-cell lineages, such as CCR7, CD28, and TNFRSF7 (CD27). A fixed-effect metaanalysis on the expression of the different modules across the datasets showed again that the 'T-cell activation' module was most significantly associated with age (Bonferroni corrected P = 3.5 9 10 À38 ) [see also Experimental Procedures]. The consistent age association of the 'T-cell activation' module, however, raises the concern that the identified modules reflect age-related changes in the proportions of cell populations in blood, as previously reported (Derhovanessian et al., 2010), rather than changes in gene expression.

Replication of coexpressed PPI modules as robust markers for aging
We conducted an independent replication study of the identified network modules as robust markers for chronological age using gene expression data from the Netherlands Twin Register and Netherlands Study of Depression and Anxiety (NTR & NESDA) consortium (N = 3535) (Boomsma et al., 2008) assayed on individuals within age range 17-79 years [Data S1, Supporting information]. An association analysis between the mean expression of a module and chronological age, adjusted for sex and the mean expression of the 'T-cell activation' module, yielded significant results for four of the five identified modules, all with directions corresponding to those found in the compendium [ Table S6, Supporting information]. These results emphasize the robustness of the findings produced by our approach and confirm that the mean module expression in whole blood of module B, C, E, and F may be considered as robust markers of chronological age.
Coexpressed PPI modules are enriched for GenAge longevity and aging genes As a validation of the identified modules, we computed whether agingrelated genes stored by GenAge (de Magalhaes & Toussaint, 2004), a database providing a comprehensive overview of aging-related genes in humans and model systems, were enriched within modules A-F (Fig. 2) [Data S1, Supporting information]. Whereas module A was supported by human derived annotations only (OR = 12.1, 95% CI 2.88-39.2, P = 6.95 9 10 À4 ), module B was solely based on knowledge derived from model organisms (OR = 16.9, 95% CI 7.26-39.1, P = 2.52 9 10 À10 ) [ Table S7, Supporting information]. Modules D, E, and F had annotations balanced over both sources, and therefore, the significance of the joint enrichment was assessed by using a resampling approach [Data S1, Supporting information], which yielded significant enrichments for modules E (P = 0.016) and F (P = 0.0029). These findings provide additional evidence that the joint expression of these modules may play a relevant role in human aging.

Module F associates with prospective survival at old age
To investigate whether the identified modules could potentially serve as biomarkers, we studied the microarray data assayed on 50 nonagenarian individuals from the Leiden Longevity Study (Passtoors et al., 2012). A left truncated Cox proportional hazard model adjusted for sex and cell counts indicates that the mean expression of module F associates with prospective survival beyond the age of 90 years (N = 50, N death = 45, HR = 0.265, 95% CI 0.12-0.57, P = 0.001). By showing that module F associates with prospective survival at old age, we illustrate its potential biological relevance. Interestingly, the ASF1A gene is part of module F and has previously been identified by our group as one of the genes that was differentially expressed in blood of members of long-lived families as compared to similarly aged controls at middle age (Passtoors et al., 2012). To confirm that the expression of the ASF1A gene in module F also associates with prospective survival at old age, we analyzed the gene expression of ASF1A measured with RT-qPCR in 74 nonagenarians from the Leiden Longevity Study (of which 24 overlapped with the micro-array experiment) for association with prospective survival. Because we observe a similar association (N death = 64, HR = 0.54, 95% CI 0.34-0.85, P = 0.008) [Fig. 3], these results indicate that modules, of which the expression in blood is consistently associated with chronological age across various datasets, may associate with variation in lifespan, and therefore provide valid gene targets for studying relevant biological endpoints in human aging.

Discussion
Age-associated changes in gene expression may provide meaningful leads to pathways affected by and involved in aging, though are generally difficult to detect consistently (de Magalhaes et al., 2009). Therefore, we constructed a large compendium of human whole blood expression studies (Goring et al., 2007;Emilsson et al., 2008;Inouye et al., 2010) comprising 2539 individuals on which we performed a novel integrative network-based analysis. This yielded fifteen consistently age-associated coexpressed PPI modules. Because the most significant age-associated module appeared to correlate with lymphocyte cell counts in an independent gene expression dataset, the expression of this module, enriched for 'T-cell activation', was subsequently used as a proxy for possible confounding shifts in the distribution of lymphocyte subsets. This enabled the identification of five age-associated modules [ Fig. 2 Panel I and IV], including three modules enriched for 'Translational elongation', 'Cytolysis' and 'DNA metabolic process' [Fig. 2B-D]. Replication in an independent cohort confirmed these findings for four of five modules [ Fig. 2B,C,E and F], underpinning the robustness of the proposed approach. The enrichments against a database for aging-related genes [Fig. 2B,E and F] emphasize the relevance of these biological findings for aging research, which is even further substantiated by the fact that the mean expression of module F associates with prospective survival at old age.

Mitochondrion-related aging
Two of the identified modules are down-regulated with age and seem to be related to the mitochondrion, though lacking any significant functional enrichment [Fig. 2E,F]. Despite the absence of functional enrichments, both modules were significantly enriched for aging-related genes, as defined by GenAge, implying that known age-related single genes can be put into a novel biological perspective by our network approach.
Module [ Fig. 2F] contains several mitochondrial factors and enzymes, like, for instance, the mitochondrial transcription termination factor MTERF, the ACADM enzyme used for fatty acid metabolism, or the mitochondrial tRNA synthetase IARS2, whose homolog was shown to increase lifespan upon disruption in worms (Smith et al., 2008). This module also includes several genes previously associated with age or age-associated diseases such as the mitotic checkpoint protein BUB3, previously associated with accelerated aging in mice (Baker et al., 2006), and the cell-cycle checkpoint protein APPBP1 found in increased quantities in the brain affected by Alzheimer's disease (Chen et al., 2003). This broad range of gene characteristics composing the module could be explained by the fact that the functionality of mitochondria is not confined to cellular energy metabolism alone, but also seems to make up an integral part of multiple cell signaling cascades including cell-cycle control and cell death (McBride et al., 2006).
Interestingly, module F also includes the ASF1A histone chaperone of which we previously have shown that its expression associates with familial longevity in the Leiden Longevity Study (Passtoors et al., 2012). We revisited the RT-qPCR data assayed on 74 nonagenarians and now show that the expression of ASF1A also associates with prospective survival. This result illustrates that modules, of which the expression in blood is consistently associated with chronological age across various datasets, may associate with variation in lifespan, and therefore provide valid gene targets for studying relevant biological endpoints in human aging.
The other mitochondrion-related module [ Fig. 2E] contains the heat shock protein HSPCA (HSP90) and the mitochondrial receptor TOMM20, which jointly play a central role in translocating preproteins into the mitochondria (Fan et al., 2011). They seem to be consistently coexpressed in blood with EIF4A2, a eukaryotic translation initiation factor Module-based meta-analysis of human aging data, E. B. van den Akker et al. 221 and DDX18, an ATP-dependent RNA helicase, of which the worm homologs were shown to extent lifespan upon disruption (Curran & Ruvkun, 2007;Smith et al., 2008). To summarize, this module seems to relate to aging by influencing protein translation and mitochondrial translocation efficiency.

Age-associated limitation of protein synthesis
One of the identified modules predominantly consisted of ribosomal proteins and translation elongation factors comprising part of the ribosomal complex [Fig. 2B]. The module was significantly enriched for 'Translational elongation' and for previous findings in model organisms with respect to aging and longevity. In addition, the module was downregulated with advancing age fitting previous observations of the aging blood transcriptome (Hong et al., 2008;Harries et al., 2011;Passtoors et al., 2012), which could be interpreted as an attempt of the cell to limit global protein synthesis in response to stress arising from damage accumulating throughout lifespan (Clemens, 2001). Whether caused by response to stress or other factors, the change in protein translation may be ascribed to the mTORC1 complex (Laplante & Sabatini, 2009). This complex modulates cellular growth and metabolisms by determining the balance between protein synthesis and degradation in response to nutrient availability. Inhibition of mTOR signaling through the mTORC1 complex not only inhibits protein synthesis, but also has been shown to positively affect the lifespan in various invertebrates and mammals (Johnson et al., 2013). Moreover, human blood transcriptome studies showed that the gene expression of mTOR pathway is down-regulated with chronological age (Harries et al., 2011;Passtoors et al., 2013) and is even associated with human familial longevity (Passtoors et al., 2012). Hence, a consistently down-regulated ribosomal module with advancing age corresponds with the age-associated demise of mTOR signaling. Although it is well established that mTOR signaling links to both lifespan regulation and 'Translational elongation', it remains to be determined whether downregulation of 'Translational elongation' is causal for human aging.

WRN-related cell-cycle checkpoint on DNA integrity
A module down-regulated with age and enriched for 'DNA metabolic process' identified in the compendium could not be replicated in the NTR & NESDA cohort [Fig. 2D]. Interestingly, this module contains the PARP1 (ADPRT) gene, which directly binds to WRN to induce apoptosis upon oxidative stress induced DNA damage and is as such a prime suspect for Werner syndrome (von Kobbe et al., 2003), a premature aging disease. Furthermore, the activity of the Parp1 protein in mononuclear cells has previously been shown to positively correlate with the species-specific lifespan across 13 mammalian species (Grube & Burkle, 1992). Taken together, findings in the compendium suggest that the lowered transcription rate of PARP1 negatively affects DNA integrity and thus lifespan, though more experiments are required to investigate this hypothesis.

Age-associated shifts in T-cell composition
Another identified module is up-regulated with age and enriched for 'Cytolysis' [Fig. 2C]. It contains several genes used to dispatch virusinfected cells and may reflect the decreased competence for fighting infections in an early stage, caused by an age-related deterioration of the immune system, known as immuno-senescence (Pawelec & Solana, 1997). We can, however, not rule out that the age-associated expression of GZMA, GZMB, and PRF1 that are part of this module point to an age-associated shift in T-cytotoxic cells (Derhovanessian et al., 2010;Napoli et al., 2012).
Though identified coregulated PPI modules may show extensive correlation with confounding factors, we should be careful to dismiss modules as such only. For instance, the 'T-cell activation' module [ Fig. 2A], which is down-regulated with age, also contained BNIP3, an inhibitor of the mTORC1 complex shown to modulate lifespan in worms, flies, and mice (Johnson et al., 2013); and FOXO1, also displaying an intricate interplay with both complexes of mTOR (Laplante & Sabatini, 2009), and shown to extent lifespan in various invertebrates (Calnan & Brunet, 2008). Additionally, human mTOR signaling may play a central role in orchestrating T-cell maturation and T-cell fate decisions (Chi, 2012), and could thereby also explain the age-associated decline in lymphocytes as marked by the 'T-cell activation' module. Taken together, these examples illustrate that what is confounding the analysis of the blood transcriptome for molecular mechanisms associated to aging is subjective to debate and might even not be possible to determine given the complex interplay between the different biological levels on which aging acts.

The proposed network approach into perspective
Network analyses have clear advantages over individual-gene analyses, as they enable the incorporation of useful prior knowledge, which can be exploited for improving the robustness of the analysis and the subsequent interpretation of the results. The improved robustness of the network approach over the individual-gene analyses was reflected by the low mutual overlap between the individual-gene results [ Fig. 1] as opposed to the high concordance between the results obtained in the compendium and replication cohort. The advantages for the interpretation were clearly illustrated by the modest insights gained from the two different strategies for individual-gene analysis ('Glycosylation site: N-linked'), as opposed to the detailed gene modules produced by our approach that can serve as a novel basis for further investigation into the molecular mechanisms underlying normative aging. Moreover, our approach is capable of inferring biological coherence from the data, without the explicit need of predefined functional groupings, as was shown by the enrichments of the identified modules found for genes within the GenAge database.
Though the analysis benefits from incorporating protein-protein interaction data, the type, and source clearly affect the results. To be as inclusive as possible for types and sources of PPI data, we have chosen to employ data obtained from the STRING database, which systematically collects and integrates interaction data derived from various sources for predicting functional relations between gene pairs. This choice results in a vast and comprehensive source of data. However, STRING data are not confined to physical interactions, as is the case with for instance IntAct (http://www.ebi.ac.uk/intact/) and unlike KEGG (http://www.genome.jp/ kegg/), STRING data are not manually curated. For network inference, a trade-off exists between the sparsity and the quality of the employed gene-gene interactions. We made use of a threshold on the quality of reported interactions that are created by STRING by benchmarking the different interaction data sources to KEGG. Varying this threshold would affect the size and nature of the obtained coexpressed PPI modules. As the threshold determines the scale of the analysis, an interesting observation is that the results can be confounded to parts of the global network that do not necessarily overlap with the predefined known biological pathways. The latter is illustrated by the fact that some of our modules are not enriched for biological pathways and could basically be valued as a strong point of our data-driven approach.

Conclusion
By applying a network approach to multiple blood transcriptomics datasets, we have identified five coexpression PPI modules that associate with chronological age in humans. The confirmation of most of our findings in an independent dataset underpins the robustness of our approach. The modules are significantly enriched for aging-related genes as curated by the GenAge database. This implies that these age-related single genes, in the absence of a clear understanding of their joint functioning belong to a network that finds its basis in protein-protein interactions and will serve as novel input for aging research. We reinforced the biological relevance of one of the modules by showing that it associates with prospective survival beyond 90 years in humans as was observed also for a single known age-related gene in this module (ASF1A). These findings collectively warrant further investigations into the biological function of module F and its potential as a biomarker for healthy aging and human longevity.

Experimental procedures
Creating the blood expression compendium Analyses were based on gene expression data derived from individuals enrolled in three large cohort studies for which details on sample inclusion and employed expression protocols are provided in depth in the original publications (Goring et al., 2007;Emilsson et al., 2008;Inouye et al., 2010). Gene expression and accompanying phenotypic data was obtained from either the original authors or from the public data repository ArrayExpress. Data quality was stringently reexamined per dataset for the presence of outlier samples or outlier measurements and annotated to a common annotation standard (EntrezGeneID). A detailed description of the data processing and an overview on the resulting sample statistics is given in the Data S1 (Supporting information) and Table 1, respectively.

Rank integration approach
A rank integration approach (Breitling et al., 2004;de Magalhaes et al., 2009) was used to identify genes consistently up-or down-regulated with age across multiple heterogeneous datasets. This type of metaanalysis integrates individual-gene statistics across datasets, by ranking the statistics per dataset and assessing the significance of the observed combined ranking using a Gamma distribution (Koziol, 2010) or through permutation. Gender adjusted linear fits between expression and age were used as gene statistics that were obtained by fitting the following multivariate linear regression model: where E ijk is the gene expression of gene i for individual j in the k th dataset, with 1 ≤ i ≤ M, 1 ≤ j ≤ N and 1 ≤ k ≤ K, where G jk and A jk are the gender and age of individual j in the k th dataset, respectively, and where e ijk is the residual error of gene i for individual j in the k th dataset. Genes were ranked on the regression coefficients between age and expression, b 2ik . The rank position of gene i in dataset k is denoted by R ik . Ranks across the datasets were integrated per gene by computing rank product statistics as previously defined by Koziol (Koziol, 2010): The significance of the observed rank products was assessed in two ways. Following Koziol, rank products RP i were transformed using: The significance of the U-statistics could be assessed by employing the gamma distribution (Koziol, 2010) or through permutation as described in the Data S1.

Extracting coexpressed PPI modules
Genes were mapped to the protein-protein interaction network (STRING v9.0, http://string-db.org/), which yielded a compendium of about 81.3% of the initial set of genes (N = 7353) in the compendium. Ranked coexpression matrices were computed for each dataset separately by computing a correlation matrix composed of first-order partial correlations between all pairs of genes adjusted for sex and subsequently assigned a rank to each of them. A higher positive correlation resulted in a higher ranking. The ranked coexpression matrices were integrated by computing rank products as in the section on individual-gene analysis. The resulting gene-gene rank product matrix together with the PPI network matrix was subsequently used as input for the method that identifies coexpressed PPI subnetworks as described in Van den Akker et al. (van den Akker et al., 2011), see also Data S1 (Supporting information). In short, a cluster analysis on the gene-gene rank product matrix yielded coexpressed modules of genes. High confidence coexpressed genes were obtained by applying a threshold on the gene-gene rank product matrix. We obtained coexpressed PPI modules by intersecting the coexpressed gene modules with the PPI network matrix. Coexpressed PPI modules were subsequently visualized using Cytoscape [Data S1, Supporting information].

Fixed-effect meta-analysis on module expressions across the blood compendium
Gene expression data were summarized per coexpressed PPI module for each dataset separately by taking the mean expression per individual over all genes in the module, resulting in a module expression for each dataset. Associations with age were tested for each coexpressed PPI module, by performing a fixed-effect meta-analysis across the four datasets using a first-order partial correlation between age and the module expression, computed with the controlling variable gender to adjust for sex differences. Per dataset k, we thus computed: q a k m k g k ¼ q a k m k À q a k g k q m k g k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À ðq a k g k Þ 2 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À ðq m k g k Þ 2 q ð4Þ where q a k m k is the correlation between age and the expression of the n th module across individuals of the k th dataset; q a k g k is the correlation between age, and gender across individuals of the k th dataset and q m k g k is the correlation between expression of the m th module and gender of individuals in the k th dataset. To correct for multiple controlling variables, higher order partial correlations were computed by repeatedly computing first-order partial correlations as described above. The function metacor of R package meta was used for integrating and testing the meta correlation statistic between age and module expression across the four datasets using default settings. Modules with significant correlations (Bonferroni corrected P-value ≤ 0.05) were considered age dependent.

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's web-site.

Fig. S1
Overview of detected consistently coexpressed PPI modules.

Table S5
Gene enrichment analyses using DAVID on the 27 consistently coexpressed PPI modules counting at least 5 genes.  Data S1 Supplemental methods.