Nội dung

Glass et al. BMC Bioinformatics (2015) 16:115 DOI 10.1186/s12859-015-0551-y RESEARCH ARTICLE Open Access A network model for angiogenesis in ovarian cancer Kimberly Glass1,2,3, John Quackenbush1,2, Dimitrios Spentzos4, Benjamin Haibe-Kains5 and Guo-Cheng Yuan1,2* Abstract Background: We recently identified two robust ovarian cancer subtypes, defined by the expression of genes involved in angiogenesis, with significant differences in clinical outcome. To identify potential regulatory mechanisms that distinguish the subtypes we applied PANDA, a method that uses an integrative approach to model information flow in gene regulatory networks. Results: We find distinct differences between networks that are active in the angiogenic and non-angiogenic subtypes, largely defined by a set of key transcription factors that, although previously reported to play a role in angiogenesis, are not strongly differentially-expressed between the subtypes. Our network analysis indicates that these factors are involved in the activation (or repression) of different genes in the two subtypes, resulting in differential expression of their network targets. Mechanisms mediating differences between subtypes include a previously unrecognized pro-angiogenic role for increased genome-wide DNA methylation and complex patterns of combinatorial regulation. Conclusions: The models we develop require a shift in our interpretation of the driving factors in biological networks away from the genes themselves and toward their interactions. The observed regulatory changes between subtypes suggest therapeutic interventions that may help in the treatment of ovarian cancer. Keywords: Network modeling, Gene regulation, Regulatory networks, Ovarian cancer, Cancer subtypes, Angiogenesis Background Ovarian cancer is the fifth leading cause of cancer death for women in the U.S. and the seventh most fatal worldwide [1]. Although ovarian cancer is notable for its initial sensitivity to platinum-based therapies, the vast majority of women eventually recur and succumb to increasingly platinum-resistant disease. Despite significant investment, improvements in patient prognosis have been slow and usually in small increments. The disease generally presents at an advanced stage (III/IV) and the five-year survival rate of advanced disease is less than 30% with median survival only slightly longer than two years [2]. Furthermore, ovarian cancer patients often undergo similar treatment regimens, mainly because the highly suspected multiple subtypes have not yet been well characterized in terms of their biological significance. * Correspondence: gcyuan@jimmy.harvard.edu 1 Dana-Farber Cancer Institute, Boston, MA, USA 2 Harvard School of Public Health, Boston, MA, USA Full list of author information is available at the end of the article In previous work, we analyzed gene expression data from 129 high-grade serous ovarian cancer samples and identified a poor-prognosis subtype characterized by the expression of angiogenic genes [3]. This subtype and the associated differences in patient survival were validated using gene expression data from a collection of 1606 ovarian cancer samples assembled from ten independent published studies. Multiple other subtypes have been proposed [4,5], but the importance of this subtyping, relative to other definitions, is that angiogenesis represents a process that is of potential clinical relevance and it is the only subtyping model which has been shown to be robust and prognostic in multiple independent datasets. Angiogenesis is one of the well-characterized hallmarks associated with cancer progression [6], playing an important role in maintaining tumor growth [7]. Angiogenesis is facilitated by interactions between cells and the extra-cellular matrix [8,9] and is associated with increased expression among a set of particular genes, © 2015 Glass et al.; licensee BioMed Central. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. 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. Glass et al. BMC Bioinformatics (2015) 16:115 including, but not limited to, matrix metalloproteinases (MMPs) [10] and VEGF [11-13]. Angiogenesis inhibition is being intensely studied as a possible therapeutic advance in ovarian cancer, but the effects in survival are still modest, suggesting that our understanding of the molecular underpinnings and biological implications of angiogenesis in this disease is still limited [2,14]. A number of clinical trials have tested drugs targeting angiogenic factors and shown these drugs have anti-cancer activity in a subset of ovarian cancer patients [15-18]; however, better understanding of the complex mechanisms driving response to these therapies is crucial to improving their efficacy and patient outcomes [17,19]. Although we have long studied ovarian cancer from the perspective of single genes and their properties, it has become clear that more integrative, systems-level analyses are necessary to better understand how the disease and its subtypes develop and progress, and how it may respond to different therapeutic interventions. The characterization of biological processes can distinguish disease states in cases where single gene biomarkers cannot [20]. The importance of applying network approaches in particular to better understand disease has previously been highlighted [21-24]. Simultaneously, it has become evident that integrative approaches that incorporate multiple sources of data to model biological systems often yield the most informative results [25-27]. Along these lines, we recently described an integrative network inference method, PANDA (Passing Attributes between Networks for Data Assimilation), that models information flow in regulatory networks by searching for agreement among various data-types, using information from each to iteratively refine predictions in the others [28]. PANDA models network interactions as communication between “transmitters” and “receivers”. In the context of PANDA’s regulatory networks, the transmitters are transcription factors and the receivers are their downstream target genes. This approach recognizes that for communication to occur, both the transmitter and the receiver have an essential role – although a transcription factor is responsible for regulating a target gene, the gene must also be available to be regulated. By constructing a “prior” regulatory network consisting of potential routes for communication (for example, by mapping transcription factor motifs to a reference genome) and integrating with other sources of regulatory information (such as protein interaction and gene expression data), one can estimate the responsibility and availability of each potential interaction, predict where communication is succeeding and failing, and deduce condition-specific network structures. Here we describe the application of PANDA to reconstructing subtype-specific regulatory networks in ovarian cancer. We identify differences in network topologies Page 2 of 17 between the angiogenic and non-angiogenic subtypes, and use this information to suggest potential therapies that may be efficacious in treating patients with angiogenicsubtype ovarian tumors. Results and discussion Building angiogenic and non-angiogenic specific regulatory network modules To begin, we downloaded mRNA expression data for 510 ovarian samples profiled by the Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp, [5]), normalized this data using fRMA [29], and mapped probes to Ensembl identifiers using the biomaRt Bioconductor package version 2.8.1 [30]. As reported in [3] we classified samples as belonging to either the angiogenic or non-angiogenic subtype; 188 were classified as part of the angiogenic subtype, and 322 were classified as nonangiogenic. We constructed separate genome-wide regulatory network models for the two subtypes. We began by mapping 111 TFs with known binding motifs to the promoters, here defined as [−750,+250] base-pairs around the transcription start site, of the 12290 genes with expression data in TCGA samples. Because transcriptional regulation involves assembly of protein complexes and allows for combinatorial regulatory processes, we collected information regarding physical protein interaction data between human transcription factors estimated using a mouse2-hybrid assay [31]. We used PANDA to integrate information from transcription factor binding motifs, protein-protein interaction data, and subtype-specific gene expression, constructing directed transcriptional regulatory networks for the angiogenic and non-angiogenic ovarian cancer subtypes (Figure 1A). For each edge that connects a transcription factor to its target gene, PANDA assigns a weight, in z-score units, that reflects the confidence level of a potential inferred regulatory relationship. Not surprisingly, we found the edge weights for the angiogenic and non-angiogenic regulatory networks to be highly correlated (Figure 1B), representing common regulatory mechanisms and processes active in both subtypes. However, we also found regulatory edges that were more strongly supported by either the angiogenic or the non-angiogenic subtype. We identified 12631 edges that we assigned to an angiogenic subnetwork, shown in red, and 15735 edges for the nonangiogenic subnetwork, shown in blue. The edges in the angiogenic subnetwork target 4081 genes while the nonangiogenic subnetwork edges target 4419; of these, 1828 genes are targeted in both subnetworks (Figure 1C), although by different upstream transcription factors. This may reflect the fact that in addition to different pathways being activated in each of the subtypes, there may also be a complex “re-wiring” of the networks around commonly targeted genes. Glass et al. BMC Bioinformatics (2015) 16:115 Page 3 of 17 Figure 1 A summary of PANDA gene regulatory network reconstruction and the identification of subtype-specific subnetworks. (A) We combine transcription factor motif and physical protein-protein interaction (PPI) data with gene expression data from each subtype to build individual network models. (B) We compare the weights of edges predicted by PANDA for each of the network models. Each point in the graph represents an edge connecting a TF to a target gene. We also define subnetworks by selecting high-probability edges specific to either the angiogenic (red) and non-angiogenic (blue) model. The number of edges and genes identified as part of these subnetworks is noted. (C) Although the subnetworks contain unique sets of edges, genes targeted by TFs in these interactions are not necessarily unique, suggesting distinct regulatory processes. Network analysis reveals biological mechanisms associated with regulatory differences We compared the corresponding subnetworks and identified a subset of transcription factors associated with the strongest regulatory changes. To do this we identified the genes targeted for regulation by a transcription factor in each of the two subnetworks and determined both an “edge enrichment” score as well as a p-value significance for the difference in the number of target genes. Specifically, we define “edge enrichment” for each transcription factor as the out-degree of that transcription factor (number of edges pointing out to a target gene) in the angiogenic subnetwork divided by the outdegree for the same transcription factor in the nonangiogenic subnetwork (multiplied by a normalization factor equal to the total number of edges in the angiogenic subnetwork divided by the total number of edges in the non-angiogenic subnetwork); the statistical significance (quantified as a p-value) of the rewiring is determined by using the hypergeometric distribution model to evaluate the overlap between edges from a transcription factor and edges specific to a particular subnetwork. On average, a given transcription factor will be associated with around 114 edges in the subtype-specific subnetworks. However some transcription factors are associated with a relatively small number (less than twenty) of edges. These transcription factors were excluded in further analysis to enhance statistical robustness. We identified ten “key” transcription factors with an edge enrichment greater than 1.5 (or less than 1/1.5) and with p < 10−3 (Figure 2A). The identified transcription factors all have established associations with angiogenesis or survival ([32-43], Figure 2B). For example NFKB1 is important for chromatin remodeling during angiogenesis [32], PRRX1 deletion causes vascular anomalies [38], and MZF1 can repress MMP2 [41], a key prognostic factor in ovarian cancer [44]. We tested if these ten transcription factors would also be identified in a simple differential expression analysis using a t-test. Only three out of ten TFs are differentially expressed at the p < 0.01 cutoff: NFKB1 (more highly expressed in the angiogenic subtype, p = 7.8 × 10−5, FDR = 4.3 × 10−4), PRRX2 (more highly expressed in the angiogenic subtype, p = 0.005, FDR = 0.015), and MZF1 (more highly expressed in the non-angiogenic subtype, p = 4.1 × 10−7, FDR = 4.2 × 10−6). Thus, PANDA identified transcription factors that are not strongly differentially expressed between the subtypes yet are known to participate in angiogenic processes. We also tested targets of each transcription factor for differential expression (see Methods) and found significant differences in target expression for six of the seven remaining TFs. For example, ARID3A is not differentially expressed between angiogenic and non-angiogenic Glass et al. BMC Bioinformatics (2015) 16:115 Page 4 of 17 Figure 2 A summary of some of the key regulatory events distinguishing the two subnetworks. (A) Transcriptional factors that have significantly more targets in one of the subnetworks compared to the other (edge enrichment > 1.5, p < 1e-3). The differential expression and methylation of each transcription factor as well as the differential expression/methylation of its target genes (in the corresponding subnetwork) is noted. Color corresponds to direction of differential expression/methylation (red: higher in angiogenic, blue: higher in non-angiogenic). (B) Genetic functions associated with these key transcription factors describing their potential role in angiogenesis and ovarian cancer. (C) A distribution of the t-statistic for differential methylation across all the genes. The shift of the distribution to the right indicates an overall increase in methylation in the angiogenic compared to the non-angiogenic samples. subtypes, but its targets have significantly (p = 2.6 × 10−22) increased expression in the angiogenic subtype. This suggests that transcription factor activity changes may not be detectable based on their own expression level increases or decreases, but that the expression of their targets can provide information on how they influence phenotype. There are many factors that could contribute to expression abnormalities in cancer, such as differences in mutations, copy-number variation and epigenetic states. We sought to determine which, if any, of these, might be contributing to the differential expression of the genes targeted by each of our ten key transcription factors. First, we investigated whether copy-number variation might explain the overall change in gene expression (see Methods). Although we find some nominally significant changes for the targets of PRRX2 (p = 0.0029) and ARID3A (p = 0.0257), these genes actually have less overall amplification in the agiogenic subtype compared to the non-angiogenic subtype (t = −2.98 and t = −2.23 respectively). This is despite the fact that their mRNA expression levels are higher in the angiogenic compared to nonangiogenic subtype. Thus copy-number variation does not appear to be the primary factor driving changes in expression of these target genes. Epigenetic factors provide an alternative mechanism for differential targeting by these transcription factors. To explore this possibility we mapped DNA methylation data from TCGA to the samples and genes used in our network reconstruction. We used a t-test to quantify any potential change in DNA methylation for each gene between the angiogenic and non-angiogenic subtypes and show the distribution of these values in Figure 2C. Overall we found DNA methylation levels to be higher in angiogenic tumor samples (mean value of t-statistic across all genes is 1.52). We next determined the differential methylation of the ten transcription factors and their targets in each of the subnetworks (Figure 2A). Compared to the rest of the genome, genes targeted by four of the ten transcription factors, ARID3A, SOX5, NKX25 and PRRX2, are associated with significantly lower Glass et al. BMC Bioinformatics (2015) 16:115 methylation in angiogenic samples. It should be noted that this lower level of methylation does not indicate hypomethylation. In fact, the average t-statistic value for the targets of these transcription factors is greater than zero (ARID3A, t = 1.17; SOX5, t = 1.31; NKX2-5, t = 0.70; PRRX2, t = 1.24), indicating that their targets, on average, have higher levels of methylation in the angiogenic compared to the non-angiogenic subtype; however, that increase in methylation is comparatively less than that experienced by genes not targeted by these regulators. We note that the TCGA methylation data was not used by PANDA to construct the networks, yet is highly concordant with the predicted patterns of gene regulation. Although many of the transcription factors that alter targeting between the aniogenic and non-angiogenic subnetworks are not significantly differentially-expressed between the subnetworks, their targets genes often are. At the same time, these target genes are also differentiallymethylated between the subtypes. Overall this analysis provides independent support of the overall network model. Both transcriptional activation and repression are used to control angiogenic pathways Transcription factors can either activate or repress gene expression. The target gene expression analysis in Figure 2A provides a preliminary indication about potential regulatory roles for the identified transcription factors. For example, although MZF1 and BRCA1 exhibit an edge-enrichment in the non-angiogenic subnetwork and are themselves also more highly expressed in the non-angiogenic samples, their targets show the opposite trend, with significantly higher expression in the angiogenic samples (p = 1.0 × 10−56 and 4.3 × 10−41, respectively). There are two scenarios consistent with these observations: (1) loss of control by these transcription factors results in the increased expression of their former targets, and (2) increased control by these transcription factors results in the decreased expression of their target genes. Interestingly, BRCA1 is known to negatively regulate IGF1 expression in breast cancer cells [43], which could inhibit angiogenesis as multiple studies have shown that increased levels of IGF1 in cancer calls leads to an increase in cell proliferation [45-47]. Similarly, MZF1 is a repressor of MMP2 [41] and is known to inhibit hematopoietic lineage differentiation in embryonic stem cells [48]. Combined with the analysis presented in Figure 2A, this suggests that MZF1 and BRCA1 act as transcriptional repressors in the nonangiogenic samples. Motivated by these observations, we classified the edges in our two subnetworks as either “activating” or “repressing” based on whether changes in the target gene's expression is correlated or anti-correlated with Page 5 of 17 subnetwork assignment. We then assigned each target gene to one of six non-overlapping classes (see Figure 3A): 1) “A+”: genes targeted only in the angiogenic subnetwork that are more highly expressed in the angiogenic subtype; 2) “A-”: genes targeted only in the angiogenic subnetwork that are more highly expressed in the non-angiogenic subtype; 3) “A+;N-”: genes targeted in both subnetworks and more highly expressed in the angiogenic subtype; 4) “N+;A-”: genes targeted in both subnetworks and more highly expressed in the non-angiogenic subtype; 5) “N-”: genes targeted only in the non-angiogenic subnetwork that are more highly expressed in the angiogenic subtype; 6) “N+”: genes targeted only in the non-angiogenic subnetwork that are more highly expressed in the non-angiogenic subtype. We used DAVID [49] to test for functional enrichment in these six classes of target genes, with the 12290 network genes taken as a background. The FDR p-values for GO “Biological Process” categories with more than one hundred members that have an FDR enrichment of less than 0.1 in at least one of our six classes of genes are illustrated as a heat map in Figure 3B. The angiogenic-activated class (“A+”) has the greatest number of significantly enriched functional categories. Many of these are associated with immune-response; processes associated with angiogenesis are also included, for example “chemotaxis,” “hematopoiesis,” “positive regulation of cell communication” and “metal ion homeostasis.” Some processes found to be enriched for genes repressed in the non-angiogenic subnetwork (“N-” genes), such as “cell adhesion” and “extracellular structure organization,” also play a role in angiogenesis. In addition, genes activated in the angiogenic subnetwork but repressed in the non-angiogenic subnetwork (“A+;N-”) include those involved in “blood vessel morphogenesis.” This suggests angiogenesis involves not only the activation of certain genes in the angiogenic subtype, but also removal of repressive regulatory interactions from the nonangiogenic subtype. In contrast, genes repressed in the angiogenic subnetwork (“A-” genes) are associated with “chromatin organization,” consistent with the observed role that epigenetics plays in distinguishing the subtypes (see Figure 2C). Genes activated in the non-angiogenic subnetwork but repressed in the angiogenic-subnetwork (“N+;A-” genes) are involved in functions such as “transcription” and “DNA metabolic process”. We also investigated whether previously identified potential biomarkers for angiogenesis were targeted in our Glass et al. BMC Bioinformatics (2015) 16:115 Page 6 of 17 Figure 3 Characteristics of six classes of differentially-targeted genes. (A) The classification of genes based on whether evidence suggests that the regulatory interactions targeting those genes are activating, repressive, or both. (B) Enriched Biological Process Gene Ontology terms (FDR < 0.1) associated with at least one of these six classes; we only included categories with at least 100 gene annotations. FDR significance is shown as a color with darker colors representing more significant enrichment. (C) Potential angiogenesis biomarkers that belong to each of the six classes of genes. Biomarkers differentially-expressed between the subtypes at an FDR < 0.1 are noted in red or blue based on whether they are more highly expressed in the angiogenic or non-angiogenic subtypes, respectively. networks, and, if so, which “class” of genes those biomarkers belonged to. In particular, we investigated thirty-five biomarkers, described in [50,51], and found twenty-two targeted in our defined subnetworks (Figure 3C). The majority (eighteen) of these biomarkers are targeted in either the “A+”, “A+;N-“or “N-“class, consistent with higher expression in the angiogenic subtype. Interestingly, many of these biomarkers are targeted in the non-angiogenic subnetwork (“A+;N-“or “N-“classes). One possible interpretation of these results is that repressive regulatory Glass et al. BMC Bioinformatics (2015) 16:115 features play a role in inhibiting angiogenic progression, in addition to transcriptional activation of these biomarkers in angiogenic tumors. Curiously, three of the four biomarkers not included in the pro-angiogenic network classes were identified in only a single study [52] that included twenty patients with inflammatory breast cancer [50]. We note that while many of the networktargeted biomarkers are also significantly differentiallyexpressed between the subtypes (FDR < 0.1 based on un-paired t-test), several are not, including VEGFA (FDR = 0.69), TP53 (FDR = 0.12), LCN2 (FDR = 0.68), KIT (FDR = 0.73) and SLC2A1 (FDR = 0.10). The identification of VEGFA and other biomarkers in our network model despite clear lack of differential-expression may indicate that our network model is able to identify important cellular regulatory alterations even in the absence of distinct changes in downstream target gene expression. Combinatorial control plays a critical role in potentiating angiogenesis Regulatory information that pertains to the core of our network can be depicted using a ring diagram representing the union of the two subnetworks (Figure 4A). In this visualization, our ten key regulators form the inner ring, while their targets, colored based on whether they Page 7 of 17 exhibit higher average expression in the angiogenic (red) or non-angiogenic (blue) samples, form the outer ring. Viewing the two subnetworks, it is clear that there is a high degree of combinatorial gene regulation. To quantify this, we applied the hypergeometric distribution model and, using the union of the genes targeted in both the angiogeneic and non-angiogenic networks as a background, tested for over-representation of genes cotargeted by specific pairs of transcription factors in the various network classifications (either “A+,” “A-,” “N+” and “N-”). Here, we focus on the three most significant pairs that include at least one of the ten identified key transcription factors. Information for all pairs can be found in the Supplemental Material (Additional file 1: Dataset S1). Note that the genes in the “A-” and “N+” classes had no combinatorial pairs significantly enriched (using a p = 10−3 cutoff). For “A+” genes, as was seen in Figure 4A, we identify significant co-regulatory associations between ARID3A, PRRX2, and SOX5 (Figure 4B-C) and these three regulators share many common targets (Figure 4D). In fact, 58% of the “A+” genes are targeted by at least one of these transcription factors, 32% by at least two, and 14% are targeted by all three, suggesting they may function as a module that coordinately regulates these genes. Figure 4 Characterizing combinatorial regulation in the subnetworks. (A) An illustration of the identified key active subnetworks. Identified key transcription factors form the inner ring and their target genes the outer ring. Target genes are colored based on whether they are more highly expressed in the angiogenic (red) or non-angiogenic subtype (blue) and are organized based on their classification. Angiogenic subnetwork edges (red) and non-angiogenic subnetwork edges (blue) extend between these rings, from the regulating transcription factor in the inner ring, to its target gene in the outer ring. (B) A table of the top three co-regulatory TF-pairs targeting “A+” genes, (C) a diagram illustrating these co-regulatory interactions, and (D) a Venn-diagram showing the overlap of the “A+” genes targeted by these TFs. (E) A table of the top three co-regulatory TF-pairs targeting “N-” genes, (F) a diagram illustrating all significant co-regulatory events between these TFs in “N-” genes, and (G) a Venn diagram showing the overlap of “N-” genes targeted by each of these TFs. Glass et al. BMC Bioinformatics (2015) 16:115 For the “N-” genes, the top three significant coregulatory transcription factor pairs include combinations of ETS1, ARNT (also known as HIF1β), MZF1 and AHR (Figure 4E). All possible pairs of these four TFs (including those that don’t include one of our “key” transcription factors) are enriched in “N-” genes at a statistically significant level (p < 1 × 10−6, Figure 4F), although MZF1 generally only shares targets of AHR or ARNT in combination with ETS1 (Figure 4G). AHR and MZF1 are among our key regulators, and, as noted previously, MZF1 is known to repress MMP2 and reduced cancer invasiveness [41]. However, ETS1 and ARNT were not among our list of key regulators, indicating that combinatorial events might be especially important for these two transcription factors. Previous reports suggest that although various ETS family members can either activate or repress angiogenic pathways [53,54], ETS1, in particular, acts as a mediator of angiogenesis [55], dimerizing with HIF2α to activate VEGFR1 and VEGFR2 [56,57]. Similarly, ARNT dimerizes with HIF1α to activate VEGF and angiogenesis [58]. However, the dimerization of AHR with ARNT inhibits ARNT/ HIF1α dimerization, thereby reducing VEGF production and subsequent angiogenesis [59]. Thus, even though ARNT/HIF1α promotes angiogenesis, the fact that ARNT/AHR dimerization inhibits angiogenesis offers an explanation for our observation that ARNT is associated with the repression of genes. Since both ETS1/HIF2α and ARNT/HIF1α interactions occur through a PAS domain [60,61], it is likely that a similar mechanism underlies our observed combinatorial enrichment of ETS1 with AHR and we hypothesize that ETS1 interaction with AHR prevents dimerization with HIFα proteins, thereby reducing VEGF production and subsequent angiogenesis. Page 8 of 17 The network model captures the effects of various treatment strategies We wished to investigate how genes identified using our network model might respond to standard or other treatment protocols. Therefore, we analyzed experimental data (GEO accession numbers GSE8057, GSE40837) measuring gene expression levels in response to several chemotherapy drugs that are commonly used to treat ovarian cancer patients and/or angiogenesis, including cisplatin, oxaliplatin and sorafenib. For each experiment, we used RMA [62] to normalize gene expression CEL files downloaded from the Gene Expression Omnibus and used a custom-CDF to map to Entrez GeneIDs [63]. We selected samples that correspond to either a treatment or control experiment and performed a t-test to quantify the differential expression of all genes between these sets of samples. Finally, we computed a summary statistic representing the aggregate differential expression value for the sets of genes within each of the six “classes” defined by our subnetworks (for more details, see Methods). The results are summarized in Figure 5A; intensity of red or blue coloration scales with the significance of increased or decreased expression, respectively, in the treatment compared to the control samples, for the genes belonging to each of our networks “classes”. Platinum-based therapies are widely used in ovarian cancer treatment regimens. Therefore, we began by investigating the effect of the chemotherapy drugs cisplatin and oxaliplatin on the expression levels of genes in A2780 human ovarian carcinoma cells (GSE8057, [64]). As a negative control, we compared expression levels of cells grown in a drug-free medium for 16 hours to their expression at 0 hours. As expected, there is little differential expression and we observe no association with any of our classes of genes (Figure 5A). We next Figure 5 Proposed therapeutic approaches. (A) A summary of the results found by comparing the expression patterns of genes in “treatment” versus “control” samples in each of the GEO datasets to “classes” of genes defined in our network analysis. We report the significance of association of differential expression with the indicated gene class, using a GSEA approach [110]; colors indicate direction of differential expression (red - increase upon treatment, blue - decrease upon treatment). (B) An illustration of some of the key findings regarding the potential mechanisms driving angiogenesis in ovarian cancer found using PANDA, as well as three potential treatments that may inhibit angiogenesis. Glass et al. BMC Bioinformatics (2015) 16:115 compared expression of cells grown for 16 hours following treatment with either cisplatin or oxaliplatin to those grown for 16 hours in drug-free medium. Curiously, in both instances we see that genes normally repressed in the non-angiogenic subnetwork (“N-”) actually increase their expression levels following treatment, suggesting that these drugs disrupt regulatory interactions that are important for repressing angiogenic activities. This is consistent with the results of a previous study [65] showing that when taken in isolation, cisplatin is not effective for treating angiogenesis in ovarian cancer. Although decreased expression in “A+” genes does not occur following platinum-based treatment, we reasoned that the effects of a VEGF-inhibiting drug should be reflected in our identified subnetworks. Biopsies of ER+ breast tumors have been collected from patients both prior to and following a clinical trial of sorafenib, and the genome-wide expression levels of genes were measured in these samples (GSE40837). Analysis of this data in the same manner as the platinum-based therapies shows a striking association with the subnetworks (Figure 5A). Genes in the “A+,” “A+;N-” and “N-” groups all show a profound decrease in expression posttreatment while genes in the “A-,” “A-;N+” and “N+” groups all increase their expression. Although this is breast rather than ovarian cancer, these results are exciting since they rely on patient samples collected from a clinical trial rather than cell-lines, illustrating a patientlevel association of an angiogenic-inhibition drug among network-defined genes. This result also serves as a positive control on our network analysis. We note that the six classes of genes we define are not wholey independent of the gene expression data, which we used both to reconstruct the networks as well as to divide target genes into distinct classes. Consequently, we would expect similar results for analyses of other “classes” of genes whose differential-expression is associated with differences between the two subtypes. However, analyzing the networks has the potential to provide additional mechanistic insight into differences between the subtypes and identify other drugs not classically associated with angiogenesis. Three treatments may synergistically inhibit angiogenic progression The optimal angiogenesis-based treatment in ovarian cancer is still a matter of ongoing investigation [14]. Commonly-used anti-angiogenic drugs mainly target VEGF, a major contributor to angiogenesis. On the other hand, as described below and illustrated in Figure 5B, several mechanisms highlighted in our network analysis suggest alternate approaches for treatment that, although speculative, could utilize existing compounds to Page 9 of 17 control, or potentially reverse, angiogenesis in ovarian cancer. For each of these proposed treatments we identified highly-related compounds and ascertained if there is a verifiable effect on gene expression in either ovarian cancer or another human system. ARNT and ETS1 dimerization with HIF1α and HIF2α, respectively, must be prevented As noted above, the dimerization of ETS1 and ARNT with HIF1α and HIF2α, respectively, generally promotes angiogenesis, although AHR may be interfering to repress target gene expression in the non-angiogenic subtype. It is therefore essential that the dimerization of these HIF proteins with ARNT and ETS1 is inhibited. HIF2α dimerizes with ARNT through a PAS-B domain, located on the C-terminus of the ARNT protein. The structural basis for this dimerization has been solved [66,67] and a small molecule ligand has been identified that dimerizes with HIF2α, decreasing affinity of the ARNT/HIF2α heterodimer [60]. Similarly, a compound has been identified that dimerizes with HIF1α, decreasing the affinity of the ARNT/HIF1α heterodimer [68]. Using either or both of these compounds, we believe one could prevent or reverse angiogenic effects driven by ARNT/HIF1α and perhaps also those driven by ETS1/HIF2α dimerization. In lieu of a dimerization inhibitor, we investigated how siRNA depletion of HIF1α affects the expression levels of genes in our identified subnetworks. We observe that in two independent experiments [69,70], “A+” genes exhibit a decrease in expression upon HIFα depletion, as we would expect from our model (Figure 5A). AHR dimerization with ARNT and ETS1 must be promoted Preventing the dimerization of ARNT/HIF1α and ETS1/ HIF2α may be insufficient; inhibition of angiogenesis is also contingent upon the dimerization of ARNT with AHR (and perhaps also ETS1/AHR). Consequently, we also suggest treatment with an AHR agonist, such as the selective AHR modular (SAhRM) 6-methyl-1,3,8-trichlorodibensofuran (6-MCDF), which has been shown to inhibit carcinogen-induced mammary tumor growth in rats [71]. TRAMP mice fed 6-MCDF in their diet had overall lower levels of serum VEFG and were five times less likely to have metastasis compared to mice on a control diet [72]. Although a known carcinogen, one of the most efficient AHR agonists is the environmental toxin 2,3,7,8Tetrachlorodibenzo-p-dioxin (TCDD). TCCD has been shown to potentiate ARNT/AHR dimerization thereby inhibiting angiogenesis and preventing vascular remodeling in rat placenta [73]. Curiously, accidental exposure to TCCD was found to potentially decrease incidence of breast and endometrial cancer in a group of women [74]. In our network, “A+” genes show a decrease in expression Glass et al. BMC Bioinformatics (2015) 16:115 upon treatment with the AHR agonist TCDD in both hepatocytes [75] and CD43+ hematopoetic cells [76] (Figure 5A). The hepatocyte data also shows a significant decrease in expression for “A+;N-” genes and “N-” genes, and an increase in expression for “A-” genes. Methylation levels across the entire angiogenic genome must be decreased In patients whose tumors have already become angiogenic, epigenetic alterations may need to be considered. One hallmark of many cancers is alteration of DNA methylation and, indeed, we found higher genome-wide methylation levels in the angiogenic subtype (Figure 2C). Interestingly, SOX5, one of our “key” transcription factors that was also found to play an important combinatorial role with ARID3A and PRRX2, contains an HMG-box which binds the minor groove [77]. Since methylation modifications occupy the major groove of DNA [78], this implies that SOX5 might be mediating key regulatory processes in the angiogenic subtype in a methylationindependent manner. This suggests that it may be necessary to decrease methylation levels across the angiogenic genome, thereby increasing competition with SOX5 binding to gene promoters, and altering their subsequent expression. This could be achieved using, for example, DNA methyltransferase (DNMT) and histone deacetylase (HDAC) inhibitors, which have already shown potential to inhibit angiogenesis in other systems [79,80]; such treatments in ovarian cancer might yield similar results. A hypomethylating agent, such as 5-azacytedine, could also be used to alter the epigenetic landscape and control angiogenic progression. With this in mind, we investigated the expression of ovarian cancer cells both prior to and post treatment with 5-azacytedine [81]; we observe a decrease in expression of the “A+” genes, consistent with our hypotheses (Figure 5A). Additional potential therapies associated with differentially-targeted genes We have identified several treatment options that target specific biological mechanisms uncovered when contrasting our network models; however, these are the result of intensive literature mining to determine suitable candidate drugs. We used the connectivity map (CMAP) [82] to determine if a gene classification based on our network analysis could also be used to identify potential drugs for treating angiogenesis in ovarian cancer. We first used genes assigned to the “A+” and “A-” classification (see Figure 3) to build a “network-signature” suitable for CMAP analysis. As expression information was included in our network model, we also built an “expression-signature” by selecting genes with the most significant changes in expression between in the angiogenic and the non-angiogenic subtypes (based on an unpaired Page 10 of 17 t-test). A cutoff of p < 1.6e-6 (FDR < 1.45 × 10−5) was used to select genes more highly expressed in the angiogenic subtype and p < 3.05e-4 (FDR < 1.4 × 10−3) to select genes more highly expressed in the non-angiogenic subtype. We used these criteria to select differentiallyexpressed genes so that the expression-signature and network-signature had equivalent dimensions (920 “up” genes and 1287 “down” genes in both cases). These two signatures share approximately 20% of their genes, with 225 in common in the “up” direction and 223 in common in the “down” direction. We used CMAP to identify drugs associated with each of these signatures. Comparing results from these two signatures will allow us to distinguish between drug candidates only identified in the network context, and those which would also be identified using a differential-expression analysis. Table 1 lists drugs significantly associated with the network-signature classification (p < 0.01). Most of these drugs are also significantly associated with the expressionsignature – not surprising given that these two signatures are non-independent; however, the ranking from the expression-signature is vastly different from that of the network-signature. One significant exception is Prestwick675, or hippeastrine, which is ranked highly in both analyses. Hippeastrine is an amaryllidaceae alkaloid with potent anti-invasive properties [83] and anticancer activities in cell lines [84]; it is also believed to contribute to the reported anti-cancer activities of the Chinese herb Lycoris aurea [85]. One of the most striking results from this analysis is that several drugs are significantly associated with the network-signature but not the expression-signature. The antitussive pentoxyverine is an agonist of the sigma-1 receptor [86], which has been shown to contribute to the induction of cancer-specific apoptosis by interleukin-24, a known inhibitor of angiogenesis [87]. Dicoumarol is an anticoagulant that is often administered to cancer patients. Anticoagulants are believed to be able to interfere with tumor angiogenesis [88] and in clinical trials their overall association with improved patient-survival, while encouraging, may to be limited to a subset of cancers [89]. Dicoumarol, in particular, inhibits furin-like activity by blocking the processing of MMP1 [34] and has been shown to abolish the TNFinduced activation of NFKB1 [90], one of our identified “core” transcription factors. Another interesting drug identified by the networksignature is harmine. Harmine was recently shown to suppress tumor growth by inhibiting angiogenic activities in endothelial cells [91] and to induce apoptosis by inhibiting the expression of MMP2 in gastric cancer [92]. Interestingly, in light of the network model we propose above, harman, a related alkaloid, stimulates AHR-dependent luciferase activity [93], and harmine is a