Network enrichment analysis: extension of gene-set enrichment analysis to gene networks

Data from a lung cancer study

For illustration we used expression data from a cohort of 123 lung cancer patients who underwent complete surgical resection at the Institut Mutualiste Montsouris (Paris, France) between 30 January 2002 and 26 June 2006. The recruitment of participants was carried out in compliance with the Helsinki Declaration, and this study was approved by the ethics committee of the Institut Gustave Roussy (Paris). All involved patients gave a written informed consent. From each patient, snap-frozen tumor and corresponding normal tissues were collected. We assayed the samples for gene expression, performed using dual-color human array from Agilent containing 41,000 gene probes; a dye-swap was employed for each sample and the log-ratio value was combined by averaging. Scanned microarray images were analyzed by using Agilent Feature Extraction software version 10.5.1.1. A loess normalization was then performed using the normalizeWithinArrays function from R package limma [16].

For a subset of 16 tumor samples we also obtained protein expression profiles. From each sample 160 μ g of protein was taken off and precipitated using four volumes of ice cold acetone. The precipitated samples were dissolved in iTRAQ dissolution buffer and digested according to manufacturers instructions (Applied Biosystems). Three proteomics platforms were used: Applied Biosystems’s MALDI TOF-TOF, Agilent’s Accurate-Mass Q-TOF and a hybrid LTQ-Orbitrap from Thermo Fischer Scientific. From these platforms, a total of 4,406 proteins were detected in at least one sample.

Altered gene sets

where log(T/N) is the log intensity-ratio of tumor vs normal expression. For the AGS with a fixed size, we considered the top 100 (k = 100) genes per individual.

The AGS from the protein expression data was defined similarly: for each sample we kept proteins whose expression exceeded the corresponding average across 16 samples. The number of proteins in the AGSs has a median of 424 (range 257 to 654).

NEA statistics and randomization

The immediate statistical question is whether an observedd _AF (k) ord _AF is significantly higher or lower than expected under the null hypothesis of no biological relation. This assessment is not trivial, as we need to respect the network topology. Crucially, it is not adequate to compare the AGS with a random list of genes. The scale-free property of most biological networks means that the degree distribution is highly skewed, with a few network hubs but a large number of sparsely connected nodes. (The degree of a node is defined as the number of neighbors connected to that node). Randomly replacing a network hub by a sparsely connected gene does not make sense. A hub gene may have a few links to genes in a FGS simply by chance, whereas for a non-hub gene the same number of links would indicate an important biological pattern.

In some simple cases, analytical calculations of the null distribution are possible, for example using the hyper-geometric distribution [17]. However, for general networks, e.g. when gene groups may share members, theoretical analyses are not feasible. In such a situation, for example, Li et al.[18] chose to limit their analysis to only cases where the gene groups did not overlap, which is a serious limitation.

A general approach in generating the null distribution of the network statistic is based on a network randomization that preserved the degree distribution. This approach was applied, for example, by [19–22]. The basic algorithm rewires the whole network by systematic swapping of network links so that the degree of each node is preserved, while its network neighbors are replaced. The algorithm is as follows:

• Step 0: Randomly select a pair of edges (A-B and C-D).

• Step 1: ‘Rewire’ the two edges in such a way that A becomes connected to D, while B connects to C. In case one or both of these new links already exist in the network, this step is aborted and a new pair of edges is selected.

Repeated applications of the above steps lead to a randomized version of the original networks. However, for a large network it is not clear at what point a sufficient level of rewiring has been achieved, and the algorithm is also inefficient as it rejects too many rewired links that are generated earlier.

The key to a fast and transparent algorithm is to first represent the network in terms of a binary adjacency table, with each node represented by a row and a column, and the cell entry ‘1’ represents connected nodes, and ‘0’ otherwise. The table is symmetric, and its margin is equal to the degree distribution of the nodes. Thus our goal is to generate a randomized table with fixed margins, similar to the hypergeometric sampling model associated with the Fisher’s exact test, but with the additional constraints of symmetry and binary entries. Adapting the hypergeometric sampling model to these constraints, the network randomization is achieved by the following sequential algorithm:

• Step 0: Start with a list of nodes 1,…,N, ordered by degreeδ₁ ≥ ⋯ ≥ δ _N > 0.

• Step 1: Assign theδ₁ links of the largest node to the other nodes randomly but with probability weighted byδ₂, …, δ _N .

• Step 2: Remove node 1, and reduce by 1 the degree of each node connected to node 1.

• Step 3: Construct a new list of remaining nodes with positive degrees, then reorder and renumber them 1,…,N, and go to Step 1. Stop if there is no remaining node with positive degree.

Likewise, the probability that the node with degreeδ _i has a link with the largest node is $\frac{{\delta }_i}{N-1}$. This is exactly the same as the weights in our algorithm. Thus, our weighting scheme simply offers a short-cut to produce truly unbiased random networks. A similar proof was used in [23]. Our proposed algorithm is computationally very efficient, since each node needs to be visited at most once. The average number of nodes that need to be visited depends on the degree distribution. In our application with ∼16,000 nodes and ∼1.5 M links (see the subsection on gene network below), this average was ∼8,200 nodes. On a 3 GhZ PC with 4 Gb RAM, using native R commands, the algorithm takes ∼60 seconds.

Denote by $n_{\mathrm{AF}}^{\ast }\left(k\right)$ the number of links between A(k) and F based on the randomized network. On replicating the randomization many times (in practice we used 100 replications), we obtain of collection of $n_{\mathrm{AF}}^{\ast }\left(k\right)$’s, which behave as a random sample from the null distribution. We can then estimateμ _AF (k) by the average of $n_{\mathrm{AF}}^{\ast }\left(k\right)$ over the randomizations, and compute $d_{\mathrm{AF}}^{\ast }\left(k\right)$ and $d_{\mathrm{AF}}^{\ast }$ based on the randomized network, and compute the right-sided p-value as the proportion of $d_{\mathrm{AF}}^{\ast }$’s ≥d _AF (similarly the left-sided p-value). The one-sided p-value is the minimum of these two, and the two-sided p-value is twice the one-sided p-value. We checked that the p-values computed for $d_{\mathrm{AF}}^{\ast }$ for a typical FGS (with median size), across all AGSs and all the randomized networks, followed a uniform distribution, indicating that the randomization worked as expected; see Figure 1.

where ${\bar{d}}_{\mathrm{AF}}\left(k\right)$ and _{s AF} (k) are the mean and standard deviation of $d_{\mathrm{AF}}^{\ast }\left(k\right)$. We use the term NEA when referring to both FNEA and MNEA. These scores can be used as a measure of activation of the FGS, thus providing a biological characterization of the AGS. Assuming a normal distribution, one could convert the scores to p-values, which are convenient for combining results of network analyses with information from other findings, such as genome-wide association studies, small-scale experiments, mutation analysis, etc.

where ${\widehat{\Pi }}_0$ is the estimated proportion of null results. Monotonicity is then imposed by taking the cumulative minimum over all p-values ≥p _k .

Comparative procedures: GEA and GSEA

where (a b c d) are the table frequencies. This can be seen as the z-statistic associated with the log odds-ratio. If some observations in the table are zero, then they are replaced with 0.5 to prevent z from being undefined.

whereg _j is the j th-ranked gene in the AGS,N _R is the number of genes in R. As with the GEA, however, GSEA also does not use any network topology in its inference. To assess the significance of observed enrichment score, a permutation test procedure was originally developed for group comparisons. In our application the AGS is defined for each person, so the permutation is performed within each person.

Gene network

We compiled a comprehensive network by merging 4 networks: (1) the data integration network of functional coupling [21], re-generated with updated datasets on protein-protein interactions (IntAct), gene expression (GEO) and sub-cellular localization (GO). This network contains 1,515,318 links, with each gene identified by the ENSEMBL gene ID; (2) curated pathways (KEGG KGML files, as of June 2010), consisting of 25,765 links; (3) known protein complexes (KEGG and CORUM, as of June 2010), consisting of 54,046 links; (4) coherent annotations (coherence was estimated by a custom software as simultaneous presence of two genes in a series of overlapping GO terms, weighted for compactness of the latter), consisting of 20,412 links.

The final network contained 1,445,027 functional links between 16,299 distinct HUPO genes (some links were lost as we translated the ENSEMBL IDs to HUPO gene symbols). Note that no evidence from the lung cancer datasets was involved in the network construction. The network is available from http://www.meb.ki.se/~yudpaw.

Functional gene sets

To characterize AGSs by involvement of known biological processes, we compiled a list of genes which are of importance in cancer:

• All 235 pathways in the KEGG database (as of 21 Apr. 2010);

• 16 GO terms that could be related to hallmarks of cancer [25];

• 7 cancer-related pathways and otherwise defined FGS from publications reporting on large-scale cancer genome projects;

• EMT, epithelial-mesenchymal transition (a pathway behind metastatic development, manual curation by S. Souchelnytskyi);

• Acidic switch (a pathway behind tumor-specific pH-shift, i.e. tumor’s ability to grow in hypoxic environment, manual curation by A. de Milito).

In total the list included 5,698 distinct HUPO gene symbols assigned to 264 FGSs (overlaps are allowed). The list of pathways with their meanings and references are given in Section B of the Additional file 1 Report. The complete list of genes for each FGS is available as an RData file at http://http:/www.meb.ki.se/~yudpaw.

Simple illustration

The principle of the network enrichment analysis is shown in Figure 2. Panel A shows a schematic diagram: in order to characterize an experimentally derived AGS with relationships to an FGS, we count the number of links between members of the two sets as given by a known background network. In this figure there are 8 links. Allowing such links expands the potential biological characterization of the altered genes using rich information on gene interaction networks. Panel B shows a realistic example from the smallest AGS from the lung-cancer proteomics data: the AGS consists of 257 deregulated proteins in one lung tumor sample (diamonds). The FGS contains 10 genes involved in epithelial-mesenchymal transition (EMT), namely TWIST1, SNAI1, SNAI2, SNAI3, SIP1, CDH1, CDH2, VIM, FN1, ACTA2. In total, 88 links were found between the AGS and FGS. Note that overlaps between AGS and FGS are allowed, and in fact these overlaps tend to produce many links; in contrast, in GEA the enrichment is based on the number of such overlaps only.

To assess the level of connectivity, we generated a series of randomized networks that preserved the degree distribution of all genes. The estimated number of links expected by chance was 51.6 and the FNEA z-score was 4.97. In the other 15 patients with proteomics data, the z-score associated with EMT ranged from 1.81 to 9.83, so that the current sample was ranked 6 out of 16. In total, 15 z-scores out of 16 exceeded 1.96. Importantly, two of the EMT proteins, VIM and FN1, were themselves found in eight and ten sample AGSs, respectively. As a result, GEA also detected some relationship between these tumors and EMT, but only 6 of the GEA scores exceeded 1.96. Looking at Figure 2B explains the higher sensitivity of NEA: network connections lead to ACTA2, CDH1, SIP1, and other EMT genes that were themselves not found in the AGS. In addition, FN1 and VIM contributed with many links.

Performance on known gene-sets

To demonstrate the biological validity of NEA, we analyzed gene sets altered in our lung cancer dataset versus known gene groups relevant to this disease. For each sample, we set the AGS to be the top 100 genes DE between tumor and normal tissues, so we had a total of 123 AGSs. Ding et al.[26] compiled a comprehensive list of 623 genes – known tumor suppressors, growth factors etc. – and screened them for somatic mutations in lung adenocarcinoma tissues. Among the 623 genes, 26 genes were declared ‘confident’ drivers in the sense that they are mutated at significantly high frequency. We performed NEA on (i) the full set of 623 genes, (ii) the set of 26 confident driver genes, and (iii) the set of 153 genes for which there were at most one mutation was detected by [26]. Another functional group was the pathway of non-small cell lung cancer from KEGG database.

Additionally, we downloaded a list of genes somatically mutated in lung cancers from the COSMIC database [27], which is potentially less biased towards previously known cancer drivers. The list of 764 genes not included in Ding et al.[26] was split into two parts: (i) 110 genes that mutated more than once and (ii) 654 genes reported only once. The latter can be regarded as a negative control. To get an unbiased comparison, we took 50 randomly sampled sets consisting of 110 genes from the 654 genes. The latter lists were expected to contain many spuriously reported passenger genes, so we expect fewer links to our lung cancer AGSs.

As for the COSMIC gene sets, even the set of multiple mutations had lower numbers of significant AGSs than the above results. Also, using FNEA, we could separate those genes that mutated once from those with more than 1 mutation, again potentially distinguishing the sets of passenger and driver mutations.

In general GEA detected only a few individuals with activated cancer pathways; see Table 2. For example, only 1 of the 123 individuals had significant enrichment of 623 cancer genes of [26], and 2 in the non-small cell lung cancer pathway. GSEA showed higher overall sensitivity compared to GEA, but GSEA appeared to be less sensitive than NEA. For example, for 26 driver genes of [26] it only identified 16 significant individuals, compared to 74 and 95 for FNEA and MNEA, respectively.

Among the most activated pathways according to GEA, we saw benzoate degradation via hydroxylation, high-mannose type N-glycan biosynthesis, biotin metabolism, lysine biosynthesis and lipoic acid metabolism. However, according to the NEA scores, few individuals were activated in these pathways. Thus, many pathways not known to be activated in cancer were found to be significant in all individuals by GEA, but not by NEA. On the other hand, putative cancer-related pathways, such as response to hypoxia, pathways in cancer, angiogenesis and the MAPK signaling pathway, were activated in most individuals according to NEA, but not according to GEA.

The results indicated that NEA can detect transcriptomic activity of cancer drivers via functional relations in the network, whereas the mutated genes themselves are not known directly. The fact that both frequently mutated genes (26 genes) and rarely mutated genes (most of the other FGSs) scored high suggests that the NEA is potentially more powerful than the mutation frequency analysis [26]. Furthermore, NEA enables investigation at the level of individuals rather than on pooled observations. Indeed, a large number of cancer drivers could not be detected by Ding et al.[26] or other studies due to low mutation frequency observed in the samples.

We further compared the procedures in terms of the false discovery rates (FDR): a better procedure should produce smaller FDR. Thus we analysed 123 individual AGSs with the 264 FGSs as described previously. Figure 3a and 3b show the average FDR curves obtained from individual AGSs. Since FNEA and GEA used a fixed number of top k genes, while MNEA and GSEA used the largest deviation over k, separate comparisons were made to compare like with like. Figure 3a indicates that for the top ranking AGS-FGS pairs, FNEA had lower FDR than GEA. Likewise, Figure 3b shows that MNEA had lower FDR than GSEA.

We also compared z-score distributions from GEA to those from MNEA and FNEA (Additional file 1: Figure S1).In general, there was little correlation between GEA and NEA scores, which would be expected as they are based on different information. In most situations NEA could achieve higher ranges of scores because there are more ways a gene can interact with other genes in a pathway than simply by being an overlap. On the other hand, the correlation between FNEA and MNEA was 0.65; since FNEA is much simpler to compute than MNEA, in practice FNEA could be sufficient for capturing the network enrichment.

Factor analysis

Once constructed, the NEA scores from each individual can be used for further analyses. For illustration we employed factor analysis [28] in order to explore the lower dimensional space that characterized molecular profiles of individual patients. The varimax rotation algorithm is used because it often provides meaningful interpretations for the factor loadings for each variable (FGS). The individual NEA scores of FGSs served as input variables in the factor analysis. In order to make the results more interpretable, the analysis was performed separately for the following three groups of pathways: (a) cancer-related pathways; (b) signaling and cancer-related pathways; and (c) metabolic pathways. These are given in Section B of the Additional file 1 Report.

First the mean and variance of FGS scores are given in Figure S2 of the Additional file 1 Report. The cancer-related pathways were neither the most enriched nor the most variable. Many other signaling KEGG pathways exceeded them in terms of these measures. The metabolic pathways seemed as variable across individuals as signaling ones. Looking at metabolic pathways most highly correlated to signaling pathways and cancer hallmarks, we found that these correlations occurred when the same enzymes were assigned to both a metabolic and a signaling pathway. For example, the RNA-specific adenosine deaminase ADAR (atrazine degradation, KEGG00710) was related to (but not a member of) the cell cycle pathway (KEGG04110) due to its molecular function. Similarly, protein phosphatases that enable biological signaling (and hence belong to multiple signaling pathways) are classified as enzymes of the inositol phosphate metabolism (KEGG00562).

Comparing the gene and protein expression data

Distinct AGSs from the transcriptomics and proteomics data allowed a systematic analysis of both similarities and differences between individual patients at molecular level. For this we compared the NEA scores of individual AGSs versus 264 FGSs for the two omics data; the plots are given in the Additional file 1: Figure S2. Both manifested higher mean scores for the five tissue growth-related pathways: adherens junction (KEGG04520), gap junction (KEGG04540), focal adhesion (KEGG04510), regulation of actin cytoskeleton (KEGG04810), and tight junction (KEGG04530). These pathways were highly enriched in the vast majority of patients. Similarly, the score distribution of metabolic pathways was quite similar between the two platforms in terms of means and variances. Intriguingly, the proteomics AGSs were strongly depleted in cancer-related FGSs, whereas in the transcription domain the cancer genesets showed the opposite pattern. We have no explanation for this observation.