Breast cancer is the malignant neoplasy with the highest incidence and mortality among women worldwide (Ferlay et al., 2014). One of the main challenges for its treatment is its heterogeneous nature, with manifestations spanning over a multitude of clinical, physiological and survival variants, resulting in differences in available therapeutic options (Polyak, 2011).

The heterogeneity of breast cancer can be traced down to the genetic level. Molecular subtyping provides a helpful tool to classify tumors by identifying common patterns in their genetic expression. Several classification algorithms that allow the classification of breast samples into molecular subtypes using different technological platforms have been developed (Perou et al., 2000; Guedj et al., 2012). A common classification scheme is given in terms of four main molecular subtypes which are luminal A, luminal B, HER2-enriched and basal-like.

The amount of genomic data available nowadays is huge and the biological information hidden in it is valuable. Using novel theoretical frameworks to analyze and understand the complex mechanisms underlying the relationships between molecules becomes a must. This is the case of the gene regulatory network approach, which uses gene expression data to model and describe the relationship between thousands of genes under particular phenotypes. In this theoretical representation, nodes (also called vertices) represent genes, while the links are an instance of interaction or relationship between those genes. These models are usually given as complex networks (Albert and Barabási, 2002; Dorogovtsev and Mendes, 2002; Newman, 2003, 2010; Boccaletti et al., 2006; Caldarelli and Vespignani, 2007).

Global organization patterns of large complex networks involve the presence of structural sub-units (subnetworks) that have been called communities or modules, broadly defined as subsets of tightly interconnected nodes so that the density of within-connections is higher than that of between-connections (Girvan and Newman, 2002; Porter et al., 2009; Fortunato, 2010). Module detection in networks, is still an open problem in computer science (Mucha et al., 2010) and there is an important variety of methods and algorithms to detect such communities (Gulbahce and Lehmann, 2008; Ahn et al., 2010; Fortunato, 2010; Xie et al., 2013; Fortunato and Hric, 2016); thus, modular structure is an issue of particular relevance, from economic and social systems (Espinal-Enríquez et al., 2015b; García-Pérez et al., 2016), to biological systems (Alcalá-Corona et al., 2016).

In the case of gene transcriptional regulatory networks (Tang et al., 2012), a module or community may correspond to a co-regulated set of genes (Wilkinson and Huberman, 2004; Zhu et al., 2008; Marbach et al., 2012; Cantini et al., 2015). In this sense, modules topology may capture some aspects of the phenomenology behind biological processes. Previous works have identified module detection as a valuable alternative for the identification of groups of genes that can associate biological features to phenotypes (Cantini et al., 2015; Alcalá-Corona et al., 2016). We followed along the lines of Alcalá-Corona et al. (2016) looking to identify biologically functional modules. We were able to find subtype-specific functional processes in modules detected by using a random-walk based community detection algorithm (Rosvall and Bergstrom, 2008).

Recently, it has been shown how the differences in transcriptional programs between breast cancer molecular subtypes are reflected in their specific transcriptional networks (de Anda-Jáuregui et al., 2016). A remaining question is whether these subtype-specific networks contain modules—i.e., subnetworks—that may be associated to known biological features. Therefore, in this work we explore the modular structure of previously inferred (de Anda-Jáuregui et al., 2016) breast cancer molecular subtype transcriptional networks for: luminal A, luminal B, basal, and HER2-enriched tumors.

For further analyses, we decided to focus in the basal subtype, since it is the one with a poorer prognosis and less options for treatment. For this subtype we identified specific communities statistically enriched for processes related to the immune system. Among these, the community containing the PSMB9 gene—composed by only 21 genes—is also the only one enriched for apoptosis-related events. Furthermore, the expression signature of molecules associated with the proteasome complex indicates that cell death-related events are strongly decreased. We suggest that communities found with this approach must be further investigated experimentally looking for alternative therapeutic options in the spirit of personalized medicine (Espinal-Enríquez et al., 2017b; Hernández-Lemus et al., 2017).

For this analysis, we used a set of 493 microarray expression profiles for breast cancer samples processed on the Affymetrix HGU133A platform. This experimental dataset was collected from the Gene Expression Omnibus (GEO) from accession numbers GSE4922 (Ivshina et al., 2006), GSE1456 (Pawitan et al., 2005), GSE7390 (Desmedt et al., 2007), GSE1561 (Farmer et al., 2005), GSE2603 (Minn et al., 2005), GSE2990 (Sotiriou et al., 2006), GSE9574 (Tripathi et al., 2008), GSE15852 (Ni et al., 2010), GSE6883 (Liu et al., 2007), and GSE3494 (Miller et al., 2005). Data preprocessing was performed following a pipeline for Robust Multi-array Average (Irizarry et al., 2003), as in Tovar et al. (2015). Breast cancer samples were classified using the well-validated PAM50 algorithm (Parker et al., 2009). PAM50 classification of this dataset was achieved previously in de Anda-Jáuregui et al. (2015, 2016).

Gene regulatory network inference from experimental data has been extensively used to unveil interactions between genes from their experimentally-measured expression levels. ARACNE (Margolin et al., 2006) is one of the most employed algorithms to calculate correlations between pairs of genes. In a nutshell, the algorithm calculates the Mutual Information (MI)—a non-parametric measure that captures non-linear dependencies between variables (Hernández-Lemus and Siqueiros-Garćıa, 2013)—in a relatively fast implementation. The method associates a MI value to each significance value (p-value) based on permutation analysis, as a function of the sample size.

In our case, we described previously the network architecture of the 4 breast cancer molecular subtypes (de Anda-Jáuregui et al., 2016). In that paper, (1) it was calculated the MI for every pair of (non-self) probe sets in the microarray platform; (2) a MI cutoff proportional to p-value, corrected for sample size was established; (3) Probe sets were mapped to HUGO gene symbols, discarding those without a corresponding gene symbol.

We also performed other three methods for network inference, based on linear correlations, namely, Pearson, Spearman and Kendall correlation measures. Additionally, we inferred mutual information-based networks by using CLR and MRNETB algorithms (Supplementary Material 1). All of these were executed with the R package minet (Meyer et al., 2008).

There are many methods and techniques for community detection (Girvan and Newman, 2002; Clauset et al., 2004; Palla et al., 2005; Newman, 2006; Rosvall and Bergstrom, 2007; Fortunato, 2010; Fortunato and Hric, 2016), we decided to use Infomap (Rosvall and Bergstrom, 2008) since it has shown to be one of the best suited algorithms for detecting modules both in performance and accuracy (Lancichinetti et al., 2009) as it was assessed in terms of the LFR benchmark (Lancichinetti et al., 2008).

In the same way in which countries, regions or cities correspond to structures with more information on a map than, say a single street; Infomap identifies the relevant structures in the network as if making a map of it from the flow of information between these structures. This is achieved by compressing the information described by a random walk inside the network, being that a random walker will spend more time in a structure with greater internal information flow in the network (a greater density of internal edges) before jumping to another, as happens when navigating between cities or countries.

The description of this walk is made using Huffman coding (Huffman et al., 1952), thus it is possible to recover the most important structures of the network by minimizing this description using information theory; just like cities that have a shorter description on a map contain much more information than a particular address whose length is longer. Hence, we applied Infomap to each of the connected components of each network, previously inferred in de Anda-Jáuregui et al. (2016), for the four breast cancer molecular subtypes—since the detection of modules only makes sense in connected networks ¹—with the aim to glimpse the modular structure of them.

Additionally, to corroborate that the network topology as well as the modular structure found with our approach are due to the intrinsic nature of tumor biology, we constructed null models for each network. In this null model, all nodes and edges for each molecular subtype network were maintained, but randomly rewired, according to the Erdös-Rényi model (Erdös and Rényi, 1959). Additionally, we constructed another null model following the method previously published in Faith et al. (2007) based on Newman et al. (2001) and Newman (2003), to preserve the same degree distribution than our inferred networks (Supplementary Material 2).

The basic hypothesis in an overrepresentation analysis (ORA), is that relevant pathways can be detected if the proportion of differentially expressed genes, within a given pathway, exceeds the proportion of genes that could be randomly expected (Garćıa-Campos et al., 2015). We performed ORA for each module, by resorting to FDR-corrected hypergeometric tests with HTSanalyzeR (Wang et al., 2011), choosing a significance below Q-value = 0.001. Categories larger than 1000 genes were discarded.

Modules were labeled according to the highest PageRank gene, i.e., each module is named by its top-PageRank centrality gene. For each subtype, we observed the number of modules and their respective size. We also observed whether or not modules are conserved across molecular subtypes to assess the existence of common regulatory process in breast cancer transcriptional programs. The hypothesis is that different network architectures determine specific community structures and hence, different processes would control the regulatory program for each molecular subtype.

We looked up for general GO categories (i.e., main branches in the ontology tree) representative of the basal subtype regulatory network. We focused on this subtype since it does not present therapeutic alternatives but cytotoxic or surgery. By identifying the processes enriched for their communities, we analyzed how specific those processes are. On the other hand, we also observed whether a process/processes appear in several modules or subtypes, providing insights regarding concomitant processes that occur in cancer, despite the network architecture (as given by the molecular subtype).

Once the functional analysis were performed we studied those genes with highest centrality measures: Betweenness centrality, clustering coefficient, shortest path length, degree and PageRank, to evaluate with more accuracy the role of those genes in the module and more importantly, in the biological process.

Finally, taking into account the differentially expressed genes, we investigated whether these expression patterns could activate or inhibit the statistically significant enriched processes. By means of z-score calculations (Krämer et al., 2014) we assessed the degree of dysregulation of said processes.

By means of mutual information-based network inference for the 4 breast cancer subtypes, we previously described (de Anda-Jáuregui et al., 2016) how differences in the transcriptional program between breast cancer molecular subtypes, reflected in different transcriptional network architectures. Table 1 briefly summarizes some of these structural differences.

From this table, it is possible to observe differences between the architectures for each network subtype. Considering these different topologies, we observe that each transcriptional network presents a characteristic modularity pattern, i.e., the gene composition of communities and the connectivity rules for each subtype are unique.

Table 1 also shows that the number of modules for each subtype is different. In Figure 2, we show a visualization of each transcriptional network, with nodes colored by module. A remaining question is whether or not those modules have a specific phenotype-dependent functional role in the regulatory program of each breast cancer molecular subtype.

With the module detection method used in this work, all network's genes are classified as an element of a module. For this reason, we have modules with associated Gene Ontology processes, but also modules which do not have statistically significant enriched categories.

Of course, different sizes of networks will have different topological parameters. In the case of mutual-information inference of networks, it is necessary to establish a threshold for valid interactions. To have large networks, a low cutoff value is needed; on the contrary, for small networks, it is mandatory to have a strict threshold of mutual information. In de Anda-Jáuregui et al. (2016) it was demonstrated that the global properties of the four networks for the breast cancer molecular subtypes are conserved throughout a large range of network sizes (3 orders of magnitude). For this reason, we performed the module detection algorithm as well as the functional enrichment with the same networks than those published in de Anda-Jáuregui et al. (2016) (top-10,000 links).

In Table 2, it can be observed how each molecular subtype has a unique set of enriched modules, regarding module sizes and number of modules, as well as the number of enriched processes for each community. It is worth mentioning that the number of enriched processes is not directly related to the size of the modules (From now on, we will name the modules according to its top-PageRank node with the nomenclature Name_of_the_Gene_comm).

For instance, LUZP4_comm for the luminal A subtype is composed by 805 genes, however, only 8 processes are enriched for this module. On the other hand, ZFP36_comm contains only 12 genes, but 12 processes are enriched on it. Processes enriched for ZFP36_comm are related to the response to stimulus and signaling.

Interestingly, each module in the luminal A network has a unique set of enriched processes. This could be related to the network structure (Figure 2A), since almost all modules are separated from each other. This is not the case for the rest of molecular subtypes. The full list of enriched processes per module for all molecular subtype networks can be found in Supplementary Material 3.

Acknowledging that the expression pattern is crucial to infer the effect of a gene regulatory network, we explored the expression profile of the COL5A2 modules in each subtype. Interestingly, for all subtypes, despite those modules contain different genes, all subtypes have an overexpression pattern, this is clearly observed in Figure 4. There, genes are colored according to their expression levels (red for overexpressed and blue for underexpressed).

As it can be observed, the COL5A2 module in each subtype reflects that those enriched processes are exacerbated, which is consistent with the fact that ECM-related processes are up-regulated, corroborating the previous observations.

We have now explored some of the differences in the regulatory programs behind common breast cancer subtypes, as seen in their underlying transcriptional networks. Based on the association between the modular structures in these networks and biological processes, we have shown that each molecular subtype has a specific functional landscape, which may be associated with the features observed in the clinical setting.

With this in mind, we will now focus on the study of the basal molecular subtype. As described in the introduction, this is the most malignant subtype of breast cancer, with the poorest prognosis, and the most restrictive therapeutic alternatives. In Figure 5A we can observe the module structure of two islands of basal subtype. Colors represent each module. We show the second largest component and the IGLC1_comm, that is not connected to the larger island, but it shares enriched processes with that island.

Figure 5B shows the enriched modules (color-filled nodes). Regarding the enriched modules found in Figure 5, we observed which GO categories are involved for each module. The results are depicted in Figure 6. There, modules are colored according to the color code of Figure 5, meanwhile labels of enriched processes are colored depending on the general type of the process.

From a visual inspection of Figure 6, it is clear that the majority of enriched categories of those modules are related to immune system processes (red labels). However, some modules have other enriched processes, such is the case of PSMB9_comm, which contains only 21 genes, but 38 GO categories are enriched by it. This is the only module in which apoptosis-related enriched processes (yellow labels) can be observed. Regarding enriched categories of PSMB9_comm we also observe:

It is also worth mentioning that 40 out of 74 enriched GO categories belong to immune system-related processes; 32 of them are unique categories for independent modules. This may imply that despite modules being sets of nodes more connected between them over the rest of the component, these communities indeed communicate with each other, and present a complete regulatory program in which specific modules act individually as a part of a whole, at least, regarding immunity. The 11 shared processes (8 of them related to immunity) at the center of the figure reinforce this suggestion.

The network modular structure found here, reflects not only how groups of genes participate in an orchestrated functional process, but also how the expression patterns of those genes are according to the specific direction in which said process functions, namely, increasing or decreasing. As an example of this last, we mention the case of PSMB9_comm in basal subtype.

As we previously mentioned, PSMB9_comm is composed by 21 genes, but that module is enriched in 38 processes, mostly related to apoptosis and immunity (Supplementary Material 7). Furthermore, network centralities in this module reveal the relevance of some genes in terms of flux of information and a coordinated regulation of sub-processes. For example, PSMB9, TAP1 and UBE2L6 are the genes with highest Betweenness centrality, clustering coefficient, degree and PageRank centralities. Because of these properties, removal of them divide the subnetwork in two modules: on the one hand, the strongly-connected module of HLA genes, mainly related to the Major Histo-Compatibility Complex (MHC), and on the other hand, the proteasome complex, mainly related to apoptosis and ubiquitination. Based on several centrality measurements, we can argue that these elements are regulating in a coordinated fashion both processes, immunity and cell death (see Supplementary Material 6).

Also important is the fact that all genes in the module are overexpressed (Supplementary Material 5). Having in mind that cell death and immunity act coordinately in basal subtype, we investigated the predicted activation or inhibition that those processes present, i.e., based on the expression profile of the geneset, which is the direction of change for a determined function. In this case, Cell death and immunity. For this purpose, we used the Diseases and Functions Analysis provided by QIAGEN's Ingenuity® Pathway Analysis (IPA®, QIAGEN Redwood City, www.qiagen.com/ingenuity), which assesses (by means of a z-score) the match between observed and predicted up/down regulation patterns, as described in Krämer et al. (2014) and Espinal-Enríquez et al. (2015a).

Interestingly, the most inhibited function predicted by the analysis is Cell death (z-scores of -8), meanwhile the most activated function, predicted by the analysis (z-scores of 7) was viral infection. This means that the expression profile of molecules related to both processes has an expression signature observed when cell death-related events are decreased and concomitantly, genes responding to a viral infection activate a defense scheme. Figure 7 shows a subset of molecules that participate in both processes and have a predicted effect on them. As it can be observed, the same expression signature, inhibits cell death and at the same time, it seems to activate viral infection genes.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.