Comparisons of gene coexpression network modules in breast cancer and ovarian cancer

Background Breast cancer and ovarian cancer are hormone driven and are known to have some predisposition genes in common such as the two well known cancer genes BRCA1 and BRCA2. The objective of this study is to compare the coexpression network modules of both cancers, so as to infer the potential cancer-related modules. Methods We applied the eigen-decomposition to the matrix that integrates the gene coexpression networks of both breast cancer and ovarian cancer. With hierarchical clustering of the related eigenvectors, we obtained the network modules of both cancers simultaneously. Enrichment analysis on Gene Ontology (GO), KEGG pathway, Disease Ontology (DO), and Gene Set Enrichment Analysis (GSEA) in the identified modules was performed. Results We identified 43 modules that are enriched by at least one of the four types of enrichments. 31, 25, and 18 modules are enriched by GO terms, KEGG pathways, and DO terms, respectively. The structure of 29 modules in both cancers is significantly different with p-values less than 0.05, of which 25 modules have larger densities in ovarian cancer. One module was found to be significantly enriched by the terms related to breast cancer from GO, KEGG and DO enrichment. One module was found to be significantly enriched by ovarian cancer related terms. Conclusion Breast cancer and ovarian cancer share some common properties on the module level. Integration of both cancers helps identifying the potential cancer associated modules. Electronic supplementary material The online version of this article (10.1186/s12918-018-0530-9) contains supplementary material, which is available to authorized users.

common. The major genes associated with susceptibility to breast and ovarian cancer are the two well-known high-penetrance cancer genes: BRCA1 and BRCA2 [1][2][3]. However, mutations in these genes account for only a very small percent of all breast cancers and ovarian cancers. Other genes such as TP53, PTEN, and STK11/LKB1, are even less common causes of breast and ovarian cancer [4]. Despite tremendous efforts to conquer such malignant diseases, research on studying the mechanism of cancer development and developing effective preventive measures is still a hot topic.
The high speed development of high-throughput technologies such as next generation sequencing of the human genome, gene expression microarrays, identification of the changes of copy numbers has dramatically accelerated the study aiming at predicting and curing such diseases. Many works have been published to address the topics on associated susceptibilities, potential biomarkers, cancer predictions and so on [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Several of them put breast cancer and ovarian cancer together in their studies [5][6][7][8][9][10]. These works either borrowed information from each other with the assumption that both cancers have similar etiologies [5,7,9], or conducted research on the differences between the pathogenic mechanisms [6,8]. Some review papers also analyzed the related research progresses of both cancers together [4,10]. In this work, we also put both cancers together to study the complex gene coexpression patterns with the network tools.
Network has been widely applied to study the complex interactions between genes, proteins, and other small molecules. It is also a popular tool for studying breast cancers and ovarian cancers [11][12][13][14][15][16][20][21][22]. When using networks to study the complex interactions, one typical concept is module, which is the densely connected subnetworks. With the network module analysis, we may infer the susceptibility genes in cancer [16], identify the biomarkers [14], and predict the prognosis of the cancer patients [11]. Current studies on cancer-related network modules are mainly based on one network, which is also true for breast cancer [14] and ovarian cancer [16]. Even in the study of more than one network, the module identification process is one network by one network, and then comparisons between different networks are performed [16]. Recently, several module identification methods on multiple networks are proposed [22][23][24]. Among them, the method proposed in [22] introduced an algorithm to find the differential modules in different networks. It mainly concentrates on the differential part. While in the paper [24], the method not only can find the modules in each network, but also can align the modules at the same time. Thus both the common and the differential parts can be detected. In the following, we compared the modules that were identified from the gene coexpression networks of breast cancer and ovarian cancer using the method in [24]. We analyzed the basic properties of the modules including density, average degree, distribution difference etc., and we did enrichment analysis of Gene Ontology (GO), KEGG pathway, Disease Ontology (DO) and Gene Set Enrichment Analysis (GSEA). By comparing the modules, both the common properties and the differences between the two cancers are detected.

Data sets
The level 3 gene expression data for breast cancer (BRCA) and ovarian cancer (OV) were downloaded from The Cancer Genome Atlas (TCGA). The gene expression data were generated with UNC AgilentG4502A. We chose the samples from the solid tissues only. There are 526 and 572 samples for both BRCA and OV, respectively. The expression value of 17,814 genes was measured. The missing data for each specific gene was imputed with the mean value of the known samples. We also downloaded the most updated protein-protein interaction (PPI) data for humans from BioGrid https://thebiogrid.org. We chose the genes having PPIs in the gene expression data, and 9603 genes were selected. To choose the differential expressed genes for BRCA and OV, t-test and Kolmogorov-Smirnov test were applied. The genes having p-values less 0.01 in both tests were chosen, where the p-values were adjusted by controlling the false discovery rate. We then got 7742 genes.

Coexpression network construction
We first computed the pairwise Pearson correlation coefficient r ij to measure the coexpression levels between gene i and gene j. To make the correlations across the two networks comparable, we 'normalized' the Pearson correlation coefficients with the same method as shown in [25]. Fisher's z transformation score for each r ij was first computed as z ij = 0.5 log 1+r ij 1−r ij . The z-scores were then normalized to make z ij follow normal distribution. After this, the 'normalized' correlations were obtained by transforming back to r ij . With these steps, the correlations in the two networks will be on the same level. We then did hard thresholding to make the network be unweighted. We took 99.5% quantile of the transformed correlation coefficients as the threshold. If the absolute value of the correlation coefficient is greater than the threshold, we assigned an edge between the corresponding genes, otherwise, there is no edge. With this method, the average degree is about 37 in both networks.

Module identification in BRCA and OV coexpression networks
Before we did module identification, we first removed the genes having no links with any other genes in both networks. The left gene number became 6779. Consider the constructed gene coexpression network G 1 (BRCA), G 2 (OV) consisting of 6779 genes. We let the adjacency matrices for both networks be A 1 , A 2 , where A k (i, j) = 1 represents there is an edge between gene i and gene j. We first applied the model proposed in [24] to cluster the genes in lower subspace for both networks. This model aims at finding the clusters in multiple networks and aligning the clusters at the same time. The main idea is to use spectral clustering to find the cluster in each network, and align them by maximizing the cluster similarity of multiple networks. Here, we only have two networks.
We assume the putative number of clusters M in both networks is given first.
Let S k be the assignment of the 6779 vertices into M clusters for the network G k , where S k im = 1, if i ∈ G k belongs to the m-th cluster, otherwise S k im = 0 for i = 1, 2, · · · , 6779, m = 1, 2, · · · , M, k = 1, 2. The optimization model is formulated as: The first term in the objective function is to do clustering in both networks separately, and the second term defines the similarity of the clusters in both networks measured by cosine function. β balances the contributions from inter-network and intra-network.
To solve the optimization problem (1), the variable S k im is relaxed. Using the same technique as spectral clustering, the optimization problem is transformed to: , and L k = D k − A k , 0 is an n × n matrix with all entries being zero. By computing the eigenvectors corresponding to the M smallest eigenvalues of matrix C, the original vertices in the networks are projected to a space of dimension M. To get the clusters, we may use k-means clustering to cluster the data points similar to spectral clustering. Due to the large size of the network, k-means does not work well. Instead, we applied hierarchical clustering with complete linkage to cluster the vertices. The distance is chosen to be the spearman distance. This is because when the size of matrix C is large, the range of the eigenvector entries is large, but their order is comparatively stable. The algorithm is summarized in 'Algorithm' .
Selection of parameter β and M The parameter β controls the connections between the vertices in both networks. When β = 0, it is equivalent to finding the clusters in two networks separately. When β becomes larger, the corresponding vertices in both networks tend to belong to the same cluster. We note that even when a group of vertices are densely connected in the first network, while their corresponding parts are isolated in the second network, the method will put all the isolated vertices in the same cluster as in the first network. Here, since both networks were controlled to have a close number of total connections, we directly set β = 1, which means the connection weight between two networks is the same as that within Algorithm: Input: Adjacency matrix A k , k = 1, 2, and M, which is the putative number of clusters.
1. Construct the matrix C ; 2. Compute the M eigenvectors u 1 , u 2 , · · · , u M corresponding to the M smallest eigenvalues of matrix C ; 3. Construct a new matrix T ∈ R 2·6779×M with columns u 1 , u 2 , · · · , u M ; 4. Cluster the points constructed from each row of matrix T into M clusters using hierarchical clustering with spearman distance.
Output: Cluster label in both networks.
both networks. We note that when β > 1 and it is within a reasonable range, the results do not change much. The number of clusters M was chosen according to the eigenvalues of matrix C. M corresponds to the first big eigengap [26]. We note that here M is not the number of clusters in either of the two networks because by choosing β, the isolated vertices in one network can also be clustered together depending on the other network. It should be the number of the union clusters in both networks. Thus the method can find both the consistent clusters and the differential clusters having quite different connection probabilities.
After the above clustering procedures, we can get a cluster label for each gene. The corresponding clusters in both networks may include different genes. We take the union of the genes as the cluster. Due to the large number of genes, the clustering method may have some bias. Some unconnected subnetworks may be clustered together. Before going to further analysis, we need to check each cluster such that it cannot include unconnected subnetworks. These resulted subnetworks are defined as the modules we identified.

Enrichment analysis
To see the associations between the identified gene modules and gene functions, we did Gene Ontology (GO) [27], KEGG pathway [28], and Disease Ontology (DO) [29] enrichment analysis for each module. GO annotates genes to biological processes (BP), molecular functions (MF), and cellular components (CC) in a directed acyclic graph structure. We only considered BP terms here. KEGG annotates genes to pathways, and Disease Ontology (DO) annotates genes with human disease associations. To see whether the identified modules show statistically significant, concordant differences between the two diseases, we also did Gene Set Enrichment Analysis (GSEA) [30]. GSEA was also done for the three types of enrichment analysis: GO, KEGG, and DO. We implemented the enrichment analysis with 'clusterProfiler' [31]. For all the cases, we let the cutoff be the Benjamini-Hochberg adjusted p-value 0.05, and recorded all the enriched terms with p-value less than 0.05.

Results
We did module identification in the gene coexpression networks of both breast cancer and ovarian cancer simultaneously. Figure 1 shows the first 350 eigenvalues of the matrix C. We chose M to be 215. After we clustered the genes into 215 clusters with hierarchical clustering, we removed those with size less than 5 and greater than 800.
With the method described, finally, we got 62 modules.
To look at the module structures in both networks, we first computed the average degree (d BRCA ,d OV ) and density(D BRCA , D OV ) for each module. Besides, we did statistical test for each module to see whether the modules in the two networks have the same connection distribution. We assume the connection between any two vertices is randomly generated following Bernoulli distribution with a given probability. We applied t-test to see whether the probability in the two networks is the same. The p-values are recorded. For all the 62 modules, after we did all the enrichment analysis, we removed those modules that have no enriched terms. Finally, 43 modules were found to have at least one type of enrichment. We put the module size, the average degree, the density, the t-test p-value, and the number of enriched terms for GO, KEGG pathway, DO, GSEA in Table 1. For GSEA, we only found GSEA GO enrichment terms.
From the t-test p-values of the modules in both coexpression networks, the structures of 29 modules are significantly different with p-values less than 0.05, of which 25 modules have larger densities in OV. Thirty one modules are enriched by GO terms, 25 modules are enriched by KEGG pathways, 18 modules are enriched by DO terms, and one module is enriched by GSEA GO terms. One module (module 14) is enriched by all the four terms. Nine modules (module 7,15,17,19,21,29,33,34,37) are enriched by GO, KEGG, and DO. We checked the details of each module and the enriched terms. All the enrichment results are put in the Additional files 1, 2, 3, 4 and 5.

GO enrichment analysis
We listed the modules that have an enriched p-value less than 10 −5 in Table 2. Three modules including module 6, 36, and 39 have the same connection distribution in both cancers. These modules have small sizes, and very significant enrichments. Six among the 7 genes in module 6 are involved in the regionalization and pattern specification process. Module 36 is mainly related to glutathione metabolic process and xenobiotic metabolic process. Five of the 6 genes are involved in these processes. Module 39 is mainly related to skeletal system development. Five among 6 genes are involved in this process. In other three modules 5, 14, and 27, the connections in the OV coexpression network are much denser. There are several isolated genes in the BRCA coexpresson network. These three modules are involved in many complex biological processes significantly. One typical example is the enriched term 'GO:0016259 selenocysteine metabolic process' . There are 89 genes in the background 16655 genes involved in this process. In module 5, 44 among the 195 (overlapping with the background) genes are involved in this process. For the term 'GO:0006614 SRP-dependent cotranslational protein targeting to membrane' , 45 genes among the 195 genes are involved in the process compared to the background 108 genes of 16655 genes in this process. These show the high correlations among the genes involved in the same process. Compared to OV, the genes involving in the same biological processes in BRCA have much less correlations, which leads to the sparser module structures.

KEGG enrichment analysis
The KEGG pathways that enrich the modules having a p-value less than 10 −4 are listed in Table 3. Module 5 is enriched by the pathway 'hsa03010 Ribosome' . Forty-six out of 130 (overlapping with the background) genes in this module belong to this pathway compared to 154 out of 7274 in the background genes. From the GO enriched terms of this module, it is clear that this module is mainly involved in the translation process. The pathways related to module 14 are mainly related to diseases. 'IL-17 signaling pathway' [32], 'TNF signaling pathway' [33], 'Rheumatoid arthritis' [34], and 'MAPK signaling pathway' [35] were shown to have relations with BRCA. For these enriched terms, no existing literatures addressed their associations with OV to the best of our knowledge. Module 23 is also related to cancers. Eight genes in this module belong to 'Pathways in cancer' . It is also enriched by 'breast cancer' with a p-value 0.049, with 3 genes associated with breast cancer in this module. This module has no enrichments related to OV. The pathway 'hsa05418 Chemical carcinogenesis' that enriches module 36 is also related to cancers. Chemical carcinogens may contribute significantly to the causation of a sizable fraction, perhaps a majority, of human cancers [36].  38. This is mainly due to the different known genes associated with the two cancers. We note that module 34 is also enriched by the female organ cancer. From all the above analysis, we found that most modules are the general modules that may not be associated with BRCA and OV. The validated enriched terms related to BRCA are much more than that of OV. One reason may be there are more researches conducted on BRCA. Among all the enriched modules, module 14 is enriched by BRCA, and module 38 is enriched by OV. Module 36 is not enriched by these two diseases, but the enriched terms are related to cancer treatment. Among these three modules, the structure of module 14 is shown to be different with a p-value 1.76E-11 in these two cancers. In the following, we give some details of these three modules.

Enrichment analysis for module 14
Module 14 consists of 25 genes. The module structure in both networks is significantly different with a p-value 1.76E-11. Figure 2 shows the module structures. It is much denser in OV compared to that in BRCA. Figure 3 shows the dotplot of GO, KEGG, DO enrichment results. For GO and DO, we plotted the first 30 enriched terms with the smallest p-values. We plotted all the 27 enriched terms for KEGG. Figure 4 shows the associations between the genes and the enriched terms. We selected the most enriched 12 GO terms with p-value less than 10 −4 , 10 KEGG terms with p-value less than 0.01, and 15 DO terms with p-value less than 0.01.
The enriched GO terms are mainly related to different responses, such at inflammatory response, response to molecule of bacterial origin, responses to wounding, immune response, etc.. Several responses have been shown to be associated with cancers [37,38]. For example, the inflammation as the seventh hallmark of cancer plays important roles in cancer development. Inflammatory cells may facilitate angiogenesis and promote the growth, invasion, and metastasis of tumor cells, which may change the genetic instability in cancer cells. Controlling the regulation of inflammatory response has a potential in both prevention and treatment of cancer [37]. There are 9 genes included in this module associated with the enriched term 'regulation of inflammatory response' , which achieves a p-value of 5.43E-08. This module also includes 7 genes that are associated with the enriched term 'adaptive immune response' , which achieves a p-value of 8.78E-05. The immune response can result in the proliferation of antigen-specific lymphocytes. When antibodies and T-cell receptors are expressed and up-regulated, immunity is acquired. Then the immune systems will initiate antigenic responses against carcinomas. A new approach to the treatment of cancer is immunotherapy, which aims to up-regulate the immune system in order that it may better control carcinogenesis [38]. Figure 4a shows the relations between the 12 GO terms with p-value less than 10 −4 and the 25 genes in this module. The genes 'CXCL1' , 'CCL20' , ' ANXA1' , 'S100A9' , 'S100A12' , 'TNF' , 'DUSP10' , 'TNFAIP2' are included in more than 6 enriched terms. There are 10 genes associated with the most enriched terms 'response to lipopolysaccharide' and 'response to molecule of bacterial origin' . For these two terms, we have not found the related literature that addresses their relations to cancers.
In the KEGG enriched terms, 'TNF signaling pathway' [33], 'Rheumatoid arthritis' [34], 'MAPK signaling pathway' [35], and 'IL-17 signaling pathway' [32] have shown to be associated with BRCA. TNF is a major inflammatory cytokine shown to be highly expressed in breast carcinomas. It induces a wide range of intracellular signal pathways including apoptosis and cell survival as well as inflammation and immunity [33]. 'MARK signaling pathway' is involved in various cellular functions, including Research on signaling pathway switch in breast cancer shows that in a large proportion of breast cancer, MARK signaling pathway is repressed, while another important pathway is activated. This mechanism may have impacts on the balance between self-renewal, proliferation, and differentiation of the tumor-initiating cells [35]. IL-17 plays crucial roles in both acute and chronic inflammatory responses. It is shown to have a direct association with breast cancer invasion in human breast tumors. IL-17 directly induced breast cancer cell invasion. There should be a potential mechanism for breast cancer invasion and tumor progression [32]. 'Rheumatoid arthritis' is mainly related to immune systems. Research shows that the risk of breast cancer is increased in non-Caucasians patients with rheumatoid arthritis while it decreased in Caucasian population [34]. Figure 4b shows the associations between the enriched KEGG terms and the genes. Four genes including 'CXCL1' , 'TNF' , 'FOS' , and 'IL1A' connect to at least 6 of the 10 terms. There are at least 6 genes associated with 'TNF signaling pathway' , 'MAPK signaling pathway' , and 'IL-17 signaling pathway' . For these enriched pathways, we have not found their associations with OV.
In the GSEA study, we ordered the genes according to the t-test p-value between the two diseases and did the analysis. Finally, 17 GO terms enrich this module. Table 5 shows the enriched terms. The 8 sequential genes having the largest t-test p-values are all in the enriched biological processes. They are 'TNFAIP3' , 'S100A9' , 'BCL3' , 'MAFF' , 'TMEM173' , 'JUNB' , 'CEBPD' , 'NFKBIZ' . Figure 5 shows the patterns for the running enrichment score. All the enriched terms have a similar pattern for these 8 genes.
By comparison of this module structure between BRCA and OV, we found module 14 is closely related to BRCA from the above enrichment analysis. From the gene coexpression network of BRCA (Fig. 2), we see some breast cancer associated genes are isolated, such as 'TNF' , 'CD55' , 'IL1A' . In OV network, these genes are highly correlated and clustered into one module. This may show that due to the tumor, the correlations of some cancer related genes decrease in BRCA.

Enrichment analysis for module 36
This module includes 6 densely connected genes in both diseases, which has a p-value 0.46 in t-test for the structure difference. Of the 6 genes, 5 genes GSTM1, GSTM2, GSTM3, GSTM4, GSTM5 encode the glutathione S-transferase that belongs to the mu class. These genes b a c Fig. 4 Enrichment results for module 14. a GO terms with p-value< 10 −4 ; b KEGG terms with p-value< 0.01; c DO terms with p-value< 0.01 function in the detoxification of electrophilic compounds such as carcinogens, therapeutic drugs, environmental toxins and products of oxidative stress, by conjugation with glutathione. Thus this module is enriched by the related biological processes such as: 'glutathione derivative metabolic process' , 'xenobiotic metabolic process' , 'sulfur compound biosynthetic process' , etc.. It is also enriched by the related pathways such as 'glutathione metabolism' , 'drug metabolism-cytochrome P450' , 'chemical carcinogenesis' , and so on. Another gene is BCAR3, which is associated with estrogen resistance and breast cancer. It is translated to the breast cancer anti-estrogen resistance protein 3. Although this module is not enriched by BRCA and OV in DO significantly, it is related to the treatment of breast cancer [39].

Enrichment analysis for module 38
Module 38 includes 11 genes. Figure 6 shows the module structure in both cancers. The connection probability in these two networks is statistically the same with t-test, although the detailed connections are different. This module is enriched by 36 GO terms with adjusted p-value less than 0.05. The gene number involved in the related biological processes is at most 4. It is enriched by 3 diseases including ovarian cancer. Four of the 14 genes are associated with OV using DO enrichment. One typical gene is 'WT1' , which connects to several other genes in the OV network, while it has no connections in the BRCA network. This gene is necessary for the development of the ovaries in females, and thus is associated with ovarian cancer. The 'WT1' protein has been found to bind a host of cellular factors such as p53. It has been ranked as the No.1 target for cancer immunotherapy by the National Cancer Institute. However, it has no associations with BRCA to the best of our knowledge. A densely connected subnetwork in this module is Kallikrein-related peptidases (KLK5, KLK6, KLK7, KLK8, KLK10). This gene family can be taken as the novel cancer biomarkers as shown in [40]. The potential of KLKs as diagnostic, prognostic, and treatment monitoring biomarkers for many types of malignancies has been extensively investigated including breast cancer and ovarian cancer. Overall, this module is closely related to cancers including breast cancer and ovarian cancer, and is more associated with ovarian cancer.

Discussion and conclusion
Breast cancer and ovarian cancer are important causes of death for women. Both of them are harmone driven and are known to have some common susceptible genes such as BRCA1 and BRCA2. Several published works have studied both cancers together with the assumption that they have the same etiologies. Coexpression network modules for both breast cancer and ovarian cancer have been studied separately by several researchers. However, there are no comparisions between the coexpression network modules between breast cancer and ovarian cancer. By comparing the modules in both cancers, we aim at finding more relations between both cancers including both the similarities and the differences.
In this work, we compared the coexpression network modules of both cancers by simultaneously identifying the modules of both cancers. By projecting the genes into a lower space representation, we did hierarchical clustering of all the genes with complete linkage of spearman distance to get the modules. 43 modules were identified to be enriched by at least one of GO, KEGG pathway, and DO terms. In most of these modules, density in OV network is larger than that in BRCA network. There are 31, 25 and 18 modules enriched by GO, KEGG pathway, and DO terms, respectively. By checking the details of each enriched term, we found that one module (module 14) is enriched by GO, KEGG and DO terms related to breast cancer. Among all GO and KEGG enriched terms, several have been validated. And there are 7 genes associated with breast cancer in this module. However, there is no enrichment information related to ovarian cancer. Another module 38, which is enriched by ovarian cancer in DO, is closely related to several cancers including breast cancer and ovarian cancer, and is more associated with ovarian cancer. From the analysis of the genes in module 36, we found that it is related to the treatment of breast cancer, although it is not enriched by breast cancer terms in DO.
Comparison of the identified modules shows the differences of breast cancer and ovarian cancer on the module level. Different modules are enriched by the two cancers significantly. It also shows some common properties of both cancers such as the KLKs family in module 36. More importantly, by simultaneous clustering of both cancers, the potential cancer related modules can be identified. For example, module 14 is significantly enriched by breast cancer related terms, but much connection information is borrowed from the ovarian cancer network. This may imply the associations between this module and breast cancer.