Global network construction.

The construction of the global network is demonstrated using gene expression, mutation, methylation, and copy number data. Input datasets are required to contain data for all samples, but can contain data different genes. Datasets can be omitted, and any data that can be converted into a binary form can be added. Indeed all continuous data are necessarily made binary (see Data section). A ‘1’ corresponds to either over- or under-expression (i.e., differential expression), copy number amplification or loss, mutation or hyper- or hypo-methylation. Conversely, a ‘0’ corresponds to normal expression, normal copy number, the absence of mutations or a normal methylation pattern. Treating over-expression of a gene identically to under-expression (and analogously, treating copy number amplification and loss, and hyper- and hypo-methylation in the same way) may seem as an oversimplification. In practice, when studying homogeneous groups of samples (e.g., belonging to an identical subtype), it is less likely that some samples in this homogeneous group would exhibit over-expression and others under-expression of the same gene. A ‘1’ in a binary dataset merely reflects an abnormal state of a gene for a particular sample, and we assume that this state of abnormality will be similar for samples exhibiting the same disease phenotype. Note that because expression data of tumor samples are analyzed relative to expression data of normal samples (see Data section), differential expression in this context does not necessarily correspond to absolute high or absolute low expression.

To construct the global network, first each binary input dataset is represented as an individual network (Fig 1a). This is achieved by converting each ‘1’ in the binary data to an undirected link connecting a sample node (labelled S_1-3 in Fig 1) with a gene node (labelled a-g). For example, if gene a is over-expressed in sample S₁, a link between node S₁ and node a is created. As a result, four separate networks (one for each input dataset) are created. Next, the individual networks are merged (Fig 1b). In the merging process, sample nodes are joined (e.g., the resulting global network only contains one sample node S₁) but gene nodes are not. Gene nodes represent an abnormal state for which the interpretation is different for each input dataset. A gene can be differentially expressed, mutated, etc. and merging genes nodes would discard this information. Consequently, the network will contain multiple gene nodes with the same gene identifier. For example, in Fig 1b, gene a is present multiple times, once as a differentially expressed gene, and once as a mutated gene.

Additionally, in order to connect heterogeneous and potentially sparse genetic aberrations with their downstream effects on gene expression, prior information in the form of known gene interactions is incorporated. The prior knowledge interaction network will ensure that a pathway is not treated as an isolated group of genes. Instead, the genes in pathway are analyzed not only relative to each other, but also relative to other pathways that may or may not overlap with the pathway under study. New nodes representing genes involved in known gene interactions are created. These ‘interaction nodes’ are connected with undirected links whenever an interaction between these nodes exists (Fig 1a and 1b). Then, the genes derived from the input datasets, representing the different abnormal states, are connected to the interaction node with the same identifier. For example, the two a nodes in Fig 1b, representing a differentially expressed and mutated state of gene a, are both connected to a single interaction node a.

Note that the newly constructed global network no longer corresponds to a physical or functional gene interaction network. Instead, the network representation is a convenient means for integrating different types of data, including prior knowledge. This is possible because all links in the network can be identically interpreted, i.e. as a qualitative “is relevant to” relation. Because the network is treated as undirected, the resulting adjacency matrix representation is symmetric. As illustrated in Fig 1c, several parts of the adjacency matrix remain empty.

Network-based similarity calculation.

The adjacency matrix representation of the global network (Fig 1c) is used to derive the similarity (a quantitative measure of relevance or importance) between samples and the gene nodes representing the abnormal states of the genes in the different datasets. The similarities are summarized in the global similarity matrix (Fig 1d). Similarity measures based on shortest path calculations would be an intuitive choice, but these measures have been shown to underperform, especially when data are qualitative, incomplete or of unknown reliability [17,18]. Because of their good performance [18–21], we used kernels calculated on graph nodes. A great multitude of kernels on graph nodes exist (see [22] for a comprehensive overview). Preliminary experiments suggested that the Laplacian Exponential Diffusion (LED) kernel yields stable results. It is calculated on the weighted Laplacian matrix L as follows [22]: Here n is the number of entities in the global network, D is the diagonal degree matrix and A(i,j) represents entry j on row i of the global network’s adjacency matrix A. K_LED(i,j) contains, at time t = α, the quantity found in node i when a unit quantity starts diffusing from node j at t = 0. The exp-operator indicates the matrix exponential. Calculating K_LED results in a similarity matrix with the same size as the original global adjacency matrix. Before continuing with the analysis, additional normalization of the similarity matrix is performed by exploiting the kernel properties of a kernel matrix K with elements k_ij:

Pathway relevance ranking.

When assessing pathway importance, the goal is to identify which pathways—containing active or inactive genes, mutated genes, genes with altered copy number or hyper- or hypo-methylated genes—are more relevant for a group of samples than others. Pathway relevance will typically correspond to abnormal behavior like activation or silencing and/or to the presence of genomic alterations. The pathways of interest should be available as sets of genes. Providing topology information for the pathways of interest is unnecessary, as the relations between genes are provided by the prior knowledge interaction network.

The pathway ranking problem can be addressed by focusing on the submatrices S_exp, S_mut, S_cnv and S_met derived from the global similarity matrix (Figs 1d and 2). These submatrices represent the similarities between the considered patient samples and genes measured in respectively the expression, mutation, copy number and methylation datasets. A high value for a sample-gene similarity indicates greater importance or relevance of a gene for that sample, and low values suggest that the gene is not important for the sample. Although each entry in these submatrices expresses a similarity between a sample and a gene present in one of the aforementioned datasets, the similarities are calculated using the entire global network. The similarity metric exploits the fact that the pathways under study are not isolated entities in the global network (note that the global network is more comprehensive than the prior knowledge interaction network). If the pathways were isolated entities, the average similarity per pathway would correspond to the average number of ‘ones’ in the pathway. Instead, pathways are interconnected through a number of mechanisms: (1) samples having ‘abnormal’ genes in different pathways or in different types of data will act as bridges between pathways and datasets, (2) pathways will typically intersect or overlap with each other and (3) the prior knowledge network also contains genes that are not present in the pathway compendium and that will act as bridges between pathways. The extra connections allow ‘abnormal’ genes from outside the pathway or from another datatype to contribute to the importance of a pathway for a particular datatype. As a result, it is possible for a mutated gene to be similar to a sample without actually being mutated in that particular sample. This can happen when the sample shares e.g. many differentially expressed genes with several other samples that do exhibit mutations in that gene.

The calculation of a single aggregate pathway score is illustrated in Fig 2. First, the matrices S_exp, S_mut, S_cnv and S_met are filtered so that they contain only genes that are present in the pathway under study. For example, gene c is not present in pathway P (containing genes a,b and d). Therefore, the third column of S_exp, corresponding to gene c, is removed. Next, for each filtered submatrix separately, the average (over all samples and genes) of the remaining similarity values in S_exp, S_mut, S_cnv and S_met is calculated, resulting in a single score per input dataset ( and ).

Because a pathway’s score will be influenced by its size, and in order to eliminate random effects, the scores-per-dataset are converted to p-values that reflect the probability of observing (at least) the obtained score purely by chance. To achieve this, a large number of random gene sets are created (in this study, we used 10,000 permutations, see implementation details below). For each random set R_i, a score is calculated for each input dataset ( and ) using the same approach as used for the pathway under study. The random score distribution obtained for each dataset then allows for the calculation of the probabilities and . Multiplying these p-values results in a single aggregate pathway score that can be used to compare pathways.

The aggregate score should not be used as an absolute measure of relevance as it does not correspond to a true joint-probability. The constituting probabilities calculated for each input dataset will very likely be correlated. Consequently, the true joint probability of observing these per-dataset similarity scores will likely be higher than the value obtained by multiplying the individual p-values. However, we can assume that the correlation between the scores obtained for each input dataset will be present for all pathways. Consequently, the aggregated score can still be used to rank pathways. It has the additional advantage that it can be broken down into four components that can be traced back to the input datasets.

Comparing the scores of pathways containing highly interconnected genes with scores obtained for random sets of, most likely, unconnected genes may appear counter-intuitive: the connectivity of the pathways might result in higher scores that could never be achieved by sets of unconnected genes. However, the p-value of both a pathway and a random set of genes is calculated using the average similarities obtained for individual genes. Those similarities in turn are a function of the abnormal states observed in the different samples (the ‘ones’ in the input datasets), and of the way the gene is connected to other relevant genes in the network. Importantly, a gene’s score is not determined by how well the gene connects to the other genes that are present in the gene set under study, be it a pathway or a random gene set. Furthermore, because all genes in the analysis are present in the prior knowledge network (genes that are not present in the prior knowledge network were filtered out), all genes have neighbors in the network. Consequently, it is perfectly possible for a gene of a random gene set to be very important, either because it contains a lot of abnormal ‘ones’, or because it lies in a network neighborhood containing other important genes that are not required to be present in the random gene set.

As an alternative to the proposed permutation strategy, one could also permute the gene labels of the input datasets, and each time measure the score for the pathway under study. However, for each permutation, the computationally very expensive calculation of the similarity matrix is required. As illustrated in S1 Fig, the probabilities obtained by the computationally expensive permutation of the dataset gene labels are highly correlated (ρ>0.99) with the probabilities obtained using the proposed method. Consequently the fast random gene set based permutation strategy was used for all subsequent analysis.

Data binarization.

In order to integrate an omics dataset in the global network representation, it needs to be converted into a list of qualitative gene-sample links or relations. Each entry in each input dataset is converted into either ‘0’ or ‘1’, representing respectively the absence or the presence of a meaningful link between a sample an a gene in the global network. Although many binarization strategies are possible, the most trivial approach was adopted: a hard threshold was set to the (absolute) values in a dataset. The dataset-dependent threshold is chosen in such a way that the fraction of entries in each data matrix that is set to ‘1’ is as a close as possible to a predefined parameter. In order to find the per-dataset threshold, a naive iterative procedure is used. This procedure is applied to all but the mutation datasets. The latter are already binary and are included without further filtering.

Gene expression, mutation, methylation and copy number data.

For breast cancer, mRNA, mutation and copy number data were downloaded from The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga) in November 2012 (corresponding to the data from [9]), and methylation data were downloaded in October 2013. 463 patients had data points for all data types and were retained in the final analysis. The ovarian cancer datasets were downloaded from TCGA in April 2014, and contained mRNA, mutation, methylation and copy number data for 447 patients. All breast and ovarian cancer data used in this study are part of the TCGA open access data tier containing only de-identified and anonymized data. For breast cancer, we focused on patients with Her2, Basal, Luminal A and Luminal B tumors. Patients with Normal-like tumors were omitted from the analysis because of the low number of available samples. All gene identifiers were mapped to Entrez gene identifiers. Genes that could not be mapped were left out of the analysis.

To reduce the problem size, all datasets were filtered to contain only genes that are present in the prior knowledge interaction network (see below). For the mRNA datasets, this resulted in a final selection of 10100 and 9463 genes for respectively breast cancer and ovarian cancer. Expression data were centered on the median of the available normal samples, and made binary according to the procedure laid out in the Method section. Rather than corresponding to absolute high or absolute low expression values, the binary data will reflect whether a gene was differentially expressed, relative to normal tissue data. No variance scaling was applied as this would increase the impact of noise for genes with very low and constant normal expression.

The mutation data were preprocessed using MutSig [13]. In order to obtain a broader selection (with a possibly high number of false positives), we used unadjusted p-values with a cut-off of 0.05. For breast cancer, the final mutation dataset contained 465 mutated genes with Entrez identifiers. For the ovarian cancer dataset, an identical procedure was employed resulting in 327 genes.

Copy number data for breast cancer and ovarian cancer were processed in the same way. Significant copy number regions were identified with GISTIC2.0 [23] using the same parameter settings as in the TCGA breast cancer overview [9]. The final dataset contained 446 and 2806 genes for respectively breast cancer, and ovarian cancer. The copy number data were made binary by applying a threshold to the absolute value of the copy number variation, according to the binarization strategy described above.

The methylation data were downloaded from TCGA and filtered using the procedure and parameters described in the ovarian cancer study of the Cancer Genome Atlas Network [16]. Summarizing, hyper-and hypo-methylated genes are kept in the data set only if the following conditions are fulfilled (parameters are for hyper-methylated genes and correspond to respectively relaxed and strict settings in the filtering procedure, for hypo-methylated genes the procedure is adjusted in a trivial way). (1) In normal tissue, the gene should only slightly be methylated (<0.5, <0.4). (2) The 90^th percentile methylation level of tumor samples should be considerably higher (0.1, 0.3) than the average normal methylation level. (3) The average gene expression of normal samples should be considerably lower (2, 3 fold) than the average gene expression of the 10% highest methylation levels in tumor samples. (4) There should exist a negative spearman correlation (<-0.2, <-0.3) between the methylation level and the gene expression level. Applying this procedure resulted in 624 hyper-methylated and 270 hypo-methylated genes for breast cancer and 42 hyper-methylated and 173 hypo-methylated genes for ovarian cancer. All methylation data are transformed from the [0,1] range to the [-0.5,0.5] range by subtracting 0.5. The methylation data were made binary according to the same binarization strategy used for the other data types.

Network data.

KEGG pathway data [24] were downloaded on November 11, 2014 using the KEGG REST-API functionality. 224 pathways (disease pathways were left out of the analysis) were merged to constitute a single network. Non-gene entities (groups and complexes) in the KEGG topology were expanded using dummy genes (with unique identifiers) that were connected to the individual genes constituting the group or complex. Interacting genes with such groups or complexes are then connected with the dummy genes rather than with the constituting genes. Because the goal is to capture and use patterns of regulation, transcription factors and their targets were explicitly added. Transcription factor—target interactions were obtained from http://encodenets.gersteinlab.org/ which is part of the ENCODE project [25]. Only interactions based on proximal TFBS data were included in the final network. The final network contained 12040 genes and 97482 interactions.

Subtype-independent pathways.

The (mitotic) cell cycle pathway is the highest scoring pathway for the Basal-like subtype, but it is also ranked high for the HER2 and the Luminal B subtypes. It is absent in the top 20 for the Luminal A subtype. Closer inspection of the components making up the aggregate score (Fig 3, S2 Fig) indicates that the score is highly determined by a substantial copy number (predominantly transcription factor MYC, cyclin CCND1 and the MDM2 oncogene) and mutation (tumor suppressor TP53) component. Mutations of TP53 are absent in Luminal A tumors, explaining the low score of this pathway for that subtype. The HER2 and the Luminal B subtype also exhibit a small methylation component that appears to be determined by the methylation status of the CCND1 and CCND2 cyclins, CDKN1C (a negative regulator of cell proliferation) and the chromatin binding MCM5 protein (S2 Fig).

The PI3K-Akt signaling pathway, responsible for cell survival and proliferation through regulation of the AKT protein kinase, scores uniformly high for all subtypes. Yet, even though the constituting score components are of a similar order of magnitude for all subtypes, the actual genetic disturbances are diverse (S3 Fig). For instance, all subtypes exhibit mutations in this pathway, but for the Basal-like tumors, this is limited to a single gene (TP53), whereas the HER2 subtype is characterized by a combination of mutations in TP53 and the PIK3CA kinase. The Luminal A and B subtypes are characterized by fewer mutations in TP53, complemented with a combination of mutations in, among others, PIK3CA, the AKT1 kinase, the PTEN phosphatase and KRAS, a kinase and well-known proto-oncogene. A similar pattern can be observed from the copy number data, where in Basal-like tumors, a uniform amplification pattern is observed for MYC and PIK3CA, whereas the other subtypes exhibit a combination of amplified MYC, CCND1, MDM2, IKBKB (a serine kinase), MCL1 (involved in apoptosis regulation) and RPS6KB1 (a protein kinase). Finally, the presence of the Basal-like—HER2 –Luminal A—Luminal B methylation gradient (S3 Fig) is reflected in the presence of a small methylation component that increases in size according to the same gradient.

Because the Jak-STAT signaling pathway, involved in cytokine and growth factor signaling, shares many genes with the PI3K-Akt signaling pathway (e.g. MYC, CCND1, MYC, AKT1, AKT2, PIK3CA, and the PIK3R1 kinase, see S4 Fig), it also scores uniformly high across all subtypes. It scores lower than the PI3K-Akt signaling pathway because TP53 is not a member of the pathway, nor many of the additional mutations (KRAS, PTEN) found in the PI3K-Akt signaling pathway. This results in a smaller mutation component in the aggregate score (Fig 1).

The TGF-Beta signaling pathway scores very high for the Basal-like and Luminal A subtypes, and intermediate (top 10) for the HER2 and Luminal B subtypes. It is involved in many aspects of development biology, and in the regulation of a wide spectrum of cellular functions, including proliferation and apoptosis [26]. Its role in breast tumor development is ambiguous as TGF-Beta signaling can both suppress and promote cancer progression [27]. Interestingly, its score is, for all subtypes, determined by a large methylation component. Indeed no other pathway has (averaged over the subtypes) so large a methylation sub-score, which is remarkable since only 6 genes of the pathway are present in the methylation dataset (S5 Fig). A possible explanation lies in the topology and size of this pathway. The TGF-Beta signaling pathway representation in the KEGG database is relatively small and linear. This makes it easy to connect upstream methylation events (almost all methylated genes in the pathway are found at the start of the signaling cascade) with downstream expression effects, resulting in large similarity scores for those genes. Together with the TGF-Beta signaling pathway, the presence of the WNT signaling pathway for all subtypes but Luminal A is notable. For the Luminal B subtype, the Hippo signaling pathway is found too. These three pathways act as large cross-talking modules [26] and it has been suggested previously that these pathways play a prominent role in (triple negative) breast cancer development [15].

The pathways described above were all ranked high for (almost) all subtypes. Functionally, these pathways take part in biological processes related to cell division, proliferation and survival. Their role in breast tumor development is well studied [9]. Interestingly, the underlying disturbances leading to the high scores can be divers. Because disturbances are connected through the a-priori network, the effect on the pathway score of a set of mutual exclusive mutations or copy number alterations occurring in that pathway can equal (and indeed exceed) the effect of a single consistent disturbance. This is clearly illustrated by the PI3K-Akt signaling pathway, where similar scores are obtained for the Basal-like subtype exhibiting only TP53 mutations and MYC amplifications, and the other subtypes where TP53 mutations and MYC amplifications are much sparser, but are complemented with other disturbances.

Subtype-specific pathways.

Because the analysis was performed separately for each subtype, the method is expected to identify not only pathways that are active or important for all subtypes, but also subtype-specific pathways. For a selection of differentially scoring pathways, Fig 4 displays a comparison of scores obtained for the different subtypes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Pathway scores compared across breast cancer subtypes for a selection of pathways.

Dark blue = Basal-like, light blue = HER2, green = Luminal A and yellow = Luminal B.

https://doi.org/10.1371/journal.pone.0133503.g004

The scores for the p53 signaling pathway are highly determined by the mutation status of TP53 (S6 Fig) and consequently, the high score for the Basal-like subtype and the low score for the Luminal A subtype are to be expected. P53 signaling is intricately related with a cell’s response to, among others, DNA damage and activated oncogenes. Differential activation of this pathway in combination with hormone receptor status is an important factor determining breast cancer progression and outcome [28]. Even though Luminal B tumors too are only infrequently mutated in TP53, their score is much higher than for Luminal A tumors. A detailed investigation of the expression data revealed that Luminal B tumors express much more genes from the p53 signaling pathway than the Luminal A tumors do (S6 Fig), suggesting that the aberrant activity of this pathway is not necessarily caused by defects in the pathway itself.

The Toll-like receptor signaling pathway, suspected to be involved in the tumor’s evasion of the immune system [29], scores very high for the Basal-like subtype, intermediate for the Luminal A subtype, and low for the HER2 and Luminal B subtypes. The differences in scores for the subtypes appear to be caused by differences in the amplification and mutation status of PIK3CA, where again the lack of PIK3CA mutations in the Basal-like subtype are compensated by additional copy number amplification of PIK3CA, copy number amplifications of the MAP3K7 and RIPK1 kinases, copy number deletions in PIK3R1, and (slight) hypo-methylation of chemokine CCL5 and toll-like receptor TLR9 (S7 Fig).

The ErbB signaling pathway couples extra-cellular growth signals to intra-cellular signaling pathways that eventually control cell survival, proliferation and motility. It is involved in tumor development and progression, and many ErbB inhibitors are currently used therapeutically [30]. Its score is highly determined by the amplification state of the ERBB2 receptor tyrosine kinase and the mutation state of PIK3CA (S8 Fig). HER2 tumors exhibit an almost total amplification of ERBB2 that translates immediately to the overexpression of ERBB2 in this subtype. Luminal A tumors lack the amplification of ERBB2 but are, like HER2 tumors, frequently mutated in PIK3CA. Luminal B tumors are infrequently mutated in PIK3CA and exhibit infrequent ERBB2 amplification, whereas in Basal-like tumors PIK3CA is not mutated, nor is ERBB2 amplified.

The MAPK-signaling pathway is one of the few pathways (together with the AMPK pathway, see below) that scores highest for the Luminal B subtype. The ample evidence that aberrant MAPK signaling promotes tumor cell proliferation, survival and metastasis, makes this pathway an interesting inhibitory drug target [31]. MAP3K1, which is known for its critical function in cell fate decisions, is one of the central genes of the pathway and is crucial in connecting upstream and downstream genes [32]. Consequently, the observed differences in scores between the subtypes are expected to depend on its mutation status. Interestingly, although MAP3K1 is more frequently mutated in Luminal A tumors (S9 Fig), the pathway score is higher for the Luminal B tumors, and the mutation sub-score for Luminal B tumors is higher than for Luminal A tumors (Fig 3). This can be due to the fact that, like for the p53 signaling pathway, the Luminal B tumors express more genes from this pathway and several genes are (slightly) more methylated (S9 Fig). The more genes that are differentially expressed or methylated, the more connections exist (through the global network) between tumor samples and mutated genes and the higher the similarity scores between the tumor samples and the mutated genes will be. Conversely, the omnipresent mutations of TP53 in the Basal-like and HER2 subtypes do not result in a particular higher score for these subtypes, since the topology of the pathway is such that the connection between TP53 and the differentially expressed and methylated genes is less straightforward [24].

Like the MAPK signaling pathway, the AMPK signaling pathway, a sensor of cellular energy status [24], scores very high for the Luminal B tumors, intermediate for Luminal A and HER2 tumors, and low for the Basal-like subtype. Here the pathway score can be explained (S10 Fig) by the inter-subtype differences in PIK3CA mutations (absent in Basal-like tumors, very frequent in Luminal A), copy number amplification of CCND1 and RPS6KB1 (unimportant in Basal-like and Luminal A, frequent in HER2 and Luminal B) and copy number loss of PIK3R1 (Basal-like).

The analysis above indicates that the presented method is able to prioritize subtype-specific pathway importance with the high score for the ErbB signaling pathway for the HER2 subtype as a typical example. The results for the p53 and MAPK-signaling pathway confirm that the integrated approach, where the analysis does not depend on a single gene nor a single type of data, is relevant and results in pathway scores that do not merely reflect mutation or copy number alteration frequencies.

Comparison with alternative approaches.

To put the proposed method into perspective, we compared it to a naïve frequency-based approach. Instead of using, for each dataset, per-gene aggregated network-based similarity scores, we used the frequency of occurrence of ‘abnormal’ states (i.e. the number of times a gene was differentially expressed, mutated, … aggregated over all samples) as the individual gene score. Using the same permutation-based approach as used for the proposed method, a p-value for each pathway can be obtained for each dataset, and the p-values can be combined to result in a single frequency-based score per pathway. This naïve approach was applied to the Luminal A breast cancer subtype, corresponding to the largest group of patients in this study. The frequency-based pathway ranking was performed on the same (filtered) datasets that were used to obtain the network-based results. The results are displayed in S2 Table. In general, the two approaches agree: even though 224 pathways were evaluated, the 50 highest scoring pathways obtained with the proposed method can almost all be found among the 50 highest scoring pathways obtained with the frequency-based approach. However, the ordering of the pathways differs considerably between the two methods. Several relevant high-scoring pathways (see Result section above) obtained with the network-based method tend to be ranked lower by the frequency-based method, as is exemplified by the MAPK signaling pathway. This clinically important pathway is ranked 14^th by the network– based method, but ranked only 32^th by the frequency based approach. The pathway exhibits, more than any other pathway, a dispersed pattern of mutual exclusive mutations. The low mutation frequencies for individual genes result in a low rank when using only frequency as a relevance criterion.

Next, we compared our pathway relevance assessment results with those mentioned in the original PARADIGM publication [6]. Re-running PARADIGM with the same data used in our experimental setup, including mutation and methylation data, was not possible. We found that the top 15 pathways mentioned in the older PARADIGM study correspond to, or overlap highly with at least one pathway of our top-ranked pathways for each subtype (S1 Table). The only exception to this is the ‘p75(NTR)-mediated signaling pathway’ ranked 6 out of 15 in the original PARADIGM study. This pathway maps to the neurotrophin signaling pathway in this study. In addition to the overlapping pathways, our approach also high-scored several other pathways (e.g. TGF-Beta, WNT and Toll-like receptor signaling). This is to be expected since firstly, in this study, additional datasets are used as input, and secondly, breast cancer subtypes are analyzed separately rather than collectively.

Important pathways common to all subtypes.

Fig 5 displays the 20 highest scoring pathways for the two most extreme subtypes. From the figure, it is obvious that in general the patient group with the worst outcome has higher pathway scores than the group with the best outcome. Only the MAPK signaling and the cell cycle pathway exhibit similar aggregate scores. The fact that many different pathways are scored differentially between the subtypes suggests that the molecular explanation for the difference in survival will most likely not have a single cause.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. The 20 highest ranking pathways for the two most extreme ovarian cancer survival-based subtypes.

The contribution of each component to the total score is indicated in a different color bar: mRNA expression (dark blue), copy number (light blue), mutation (green) and methylation (yellow).

https://doi.org/10.1371/journal.pone.0133503.g005

In previous studies, as summarized in [16], a number of important deregulated pathways have been related to ovarian cancer: RB-signaling (cell cycle), RAS/PI3K signaling (PI3K-Akt signaling), FOXM1 (not present in the KEGG database, but overlapping with the cell cycle pathway) and Notch signaling. In this study, the highest scoring pathways are the MAPK signaling pathway, the cell cycle pathway, the PI3K-Akt signaling pathway and the WNT signaling pathway. The absence of the Notch signaling pathway in the top 20 of highest scoring pathways for either subtype is striking. A closer inspection of the actual data (S11 Fig) reveals that the Notch signaling pathway harbors few (epi-) genetic disturbances, both in terms of disturbances per gene as well as in terms of the number of disturbed genes. Additionally, as the KEGG representation of the Notch signaling pathway only contains 24 genes, it is possible that the pathway information is incomplete. However, a close cooperation between MAPK signaling and Notch signaling has been described in other tumors [33], and MAPK signaling is top-scored for all subtypes. The presence of the WNT signaling pathway, implicated previously in ovarian cancer tumorigenesis [34,35], is noteworthy. Its high score will likely be determined (S12 Fig) by TP53 mutations and copy number amplification in e.g. MYC, the NDK2 kinase and WNT5B (encoding for a signaling protein).

Differential subtype analysis.

More than in the pathways that score uniformly high across subtypes, we are interested in pathways that score differently for the survival-based patient groups. To achieve that, we focused on the 20 highest scoring pathways for the bad-outcome group. For each of these pathways, we calculated the ratio of the score for the bad-outcome group and the corresponding score in the good-outcome group (Fig 6). The assumption is that the larger a score ratio deviates from one, the more likely a pathway is involved in a process that determines outcome.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Ratio of bad-outcome pathway scores and the corresponding good-outcome scores.

A ratio of ‘1’ indicates that the pathway scores equally high for patients in the bad-outcome group and patients in the good-outcome group. Values larger than 1 indicate higher pathway importance / activity for the bad-outcome group. Pathways shown are limited to the top-20 highest scoring pathways in the bad-outcome group.

https://doi.org/10.1371/journal.pone.0133503.g006

The scores of the pathways in the bad-outcome group are consistently higher than the corresponding scores in the good-outcome group, with the exceptions of the MAPK signaling and cell cycle pathways yielding score ratios respectively equal to 0.97 and 0.94. The highest score ratios are found for Jak-STAT signaling, ECM receptor interaction, cAMP signaling, TGF-Beta signaling and estrogen signaling. S14–S18 Figs illustrate that single (epi-)genetic aberrations cannot be easily linked to the differential scoring of these 5 pathways, except for slight increases in MYC, CACNA1D (a Calcium channel) and KRAS amplification, and infrequent CREBBP (a transcription factor activator) mutations. Still, all five pathways are involved in tumor development and proliferation. Jak-STAT signaling is since long linked to cell migration and survival in ovarian cancer [36] and often, the expression of, or the responsiveness to TGF-Beta is lost, leading to cell proliferation [37]. Furthermore, it is known that the disruption of the extra cellular matrix (ECM interaction) can lead to the overproduction of growth factors that promote abnormal proliferation [38]. Estrogen signaling is thought to be involved in the establishment of a pro-tumorigenic micro-environment, and the long-term use of estrogen-only hormone replacements is linked to ovarian cancer development [39]. Finally, cAMP signaling is involved in cell survival in ovarian cancer and the presence of mutant CREBBP proteins can lead to tumorigenesis [40].

Interestingly, when all pathways are considered instead of only the pathways that are important for the bad-outcome patient group, the ratio score ranking (S13 Fig) is topped by the Notch signaling pathway event though this pathway is not high-scored for any patient group. As mentioned above, Notch signaling is involved in embryonic development, proliferation, differentiation and apoptosis, and is implicated in ovarian (and many other) cancers [16,41,42]. Its large score ratio could be attributed to the differential amplification of DLL3 (S11 Fig), which recently has been reported to be the second most altered gene in the NOTCH3 pathway [43].