Genome-wide identification of key modulators of gene-gene interaction networks in breast cancer
- Yu-Chiao Chiu†1, 2,
- Li-Ju Wang†3,
- Tzu-Hung Hsiao4Email author,
- Eric Y. Chuang2, 5Email author and
- Yidong Chen1, 6Email author
© The Author(s). 2017
Published: 3 October 2017
With the advances in high-throughput gene profiling technologies, a large volume of gene interaction maps has been constructed. A higher-level layer of gene-gene interaction, namely modulate gene interaction, is composed of gene pairs of which interaction strengths are modulated by (i.e., dependent on) the expression level of a key modulator gene. Systematic investigations into the modulation by estrogen receptor (ER), the best-known modulator gene, have revealed the functional and prognostic significance in breast cancer. However, a genome-wide identification of key modulator genes that may further unveil the landscape of modulated gene interaction is still lacking.
We proposed a systematic workflow to screen for key modulators based on genome-wide gene expression profiles. We designed four modularity parameters to measure the ability of a putative modulator to perturb gene interaction networks. Applying the method to a dataset of 286 breast tumors, we comprehensively characterized the modularity parameters and identified a total of 973 key modulator genes. The modularity of these modulators was verified in three independent breast cancer datasets. ESR1, the encoding gene of ER, appeared in the list, and abundant novel modulators were illuminated. For instance, a prognostic predictor of breast cancer, SFRP1, was found the second modulator. Functional annotation analysis of the 973 modulators revealed involvements in ER-related cellular processes as well as immune- and tumor-associated functions.
Here we present, as far as we know, the first comprehensive analysis of key modulator genes on a genome-wide scale. The validity of filtering parameters as well as the conservativity of modulators among cohorts were corroborated. Our data bring new insights into the modulated layer of gene-gene interaction and provide candidates for further biological investigations.
As technologies of high-throughput profiling advance, a large volume of post-transcriptional gene interaction maps has been established. For instance, the Kyoto Encyclopedia of Genes and Genomes (KEGG) is a knowledge-based curation of abundant genomic pathways among species . Such maps provide better understanding to the molecular signaling in cells, however, they are typically derived under a certain cellular condition in a single model cell line. In light of the dynamicity and complexity of gene interactions (reviewed in [2, 3]), a higher-order layer of interaction networks that considers gene-gene relationships modulated by (i.e., dependent on) key modulator genes, namely modulated gene interaction, was proposed (reviewed in ). In this sense, interaction of two genes can be strengthened specifically when a modulator gene is expressed at high or low abundance. The scenario provides flexibility and interpretability to condition-specific and dynamic interaction networks.
In breast cancer, estrogen receptor (ER) is the best-studied modulator gene. It governs the coexpression among several keratin genes in breast cancer patients . Also, topological and temporal changes were observed in a transcription factor interaction network of MCF7 cells upon 17β-estradiol stimulation . A comprehensive in silico investigation revealed compact gene-gene and function-function interaction networks modulated by ER and discovered the prognostic value of ER-modulated interaction between TGFβ and NFκB . By a co-modulation analysis, we previously showed ten experimentally chosen genes jointly modulated up to two-thirds of all gene pairs, with an implication in cellular processes associated with hormone stimulus . Taken together, these reports demonstrate the existence and functional significance of modulated gene interaction, and motivate a comprehensive search for key modulator genes. Based on mutual information, a modulator inference by network dynamics (MINDy) was developed to systematically identify modulators of transcription factor (TF)-target gene interactions . However, due to a heavy computational burden caused by permutation-based assessment of statistical significance, the method was limited to the investigation of specific TFs and relied on prior knowledge of TF-target relationships. Recently, we exploited the transformability of Pearson correlation coefficients to devise a highly efficient modulated gene/gene set interaction (MAGIC) analysis and realized the exploration into genome-wide interaction networks modulated by a modulator gene . However, a reverse-engineering study for a genome-wide identification of key modulators is still lacking.
In the present study we proposed a systematic workflow that incorporates the MAGIC algorithm to analyze gene expression profiles of breast tumors. Comparing samples with high and low expression levels of a modulator, four modularity parameters were designed to measure modulator-dependent changes in gene interaction at two layers. One was focused on the summary of genome-wide changes, while the other assessed the scale and information flow in the core subset of modulated interaction pairs. Genes with significantly high values of parameters were defined as key modulators and validated by three independent cohorts. Functional annotation analysis was performed to study the functional involvement of these modulators. Collectively, this report describes a novel genome-wide search for key modulators in breast cancer and unveils the functional landscape of modulated gene interactions.
We downloaded and reanalyzed four public gene expression microarray datasets of breast cancer patients from the Gene Expression Omnibus database  and The Cancer Genome Atlas (TCGA). A dataset of 286 lymph-node negative breast tumors (GSE2034) , profiled by Affymetrix Human Genome U133A Arrays, was analyzed for the identification of key modulator genes. We validated the findings in three large independent cohorts, GSE2990 , GSE4922 , and TCGA [14, 15]. Gene-level intensity values of GSE2034, GSE2990, and GSE4922 were calculated by reprocessing of Affymetrix CEL files by Robust Microarray Analysis (RMA) algorithm, representation of each gene by the most informative probe (measured by the coefficient of variation (CV)), and removing non-informative genes, as previously described [7, 8]. For the TCGA dataset, we used pre-normalized level-3 (gene-level) data.
Analysis of overall changes in genome-wide gene interactions networks
Analysis of core modulated gene interactions networks
where N ′ denotes the assigned sample size and was set equal to N/4 in this study for the two groups were equally sized.
By the two MAGIC parameters, core m-modulated pairs were extracted and merged into a m-modulated gene interaction network. Here we defined three more modularity parameters to measure the size and information flow of the network, Parameter 2: numbers of nodes (genes), Parameter 3: numbers of edges (gene interactions), and Parameter 4: connectivity (average node degree). Statistical significance of the three parameters were tested by 10,000-time random permutations as described in Eq. 5.
Functional annotation analysis and network visualization
Functional annotation analysis was performed by the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 [16, 17] to analyze the enrichment of key modulator genes in biological functions and processes. We focused on Gene Ontology (GO) terms of molecular functions, biological processes, and cellular components. We used the Functional Annotation Clustering tool to group GO terms to eliminate potential biases from highly similar terms. Gene interaction networks were analyzed and visualized by an open source software Cytoscape v3.2.1 , with nodes and edges representing genes and gene interactions, respectively, and node size denoting node degree.
Results and Discussion
Genome-wide identification of key modulator genes
The present study is aimed to systematically screen for key modulator genes from global gene expression data. As illustrated in Fig. 1, we selected and compared the samples with high (top 25%) and low (bottom 25%) expression of a candidate modulator gene m. Four parameters were designed to measure the modularity of m from two aspects, one at a genome-wide level and the others focusing on the core subnetwork only. ACI score (parameter 1) represents the overall change in interaction strengths between genome-wide gene interaction networks formed in the two sample groups. Focusing on the core sub-network (m-modulated gene interaction network) constructed merely by significantly changed edges, we further designed three parameters (namely, number of nodes, number of edges, and connectivity) to quantify the scale and information flow affected by the modulation of m. For each m, significance of the four parameters was tested by random permutation of dataset. Mathematical details are described in the Methods section.
Properties of modularity parameters
Identification and validation of key modulator genes
Top 20 modulator genes
SFRP1-modulated gene interaction network
Top 6 clusters of GO terms enriched in SFRP1-modulated gene interaction network
Cluster 1 (enrichment score: 10.71)
proteinaceous extracellular matrix
extracellular matrix part
extracellular matrix structural constituent
Cluster 2 (enrichment score: 7.13)
response to organic substance
response to hormone stimulus
response to steroid hormone stimulus
response to endogenous stimulus
response to estrogen stimulus
Cluster 3 (enrichment score: 6.56)
extracellular region part
Cluster 4 (enrichment score: 6.27)
cell cycle process
cell cycle phase
mitotic cell cycle
Cluster 5 (enrichment score: 5.80)
microtubule cytoskeleton organization
Cluster 6 (enrichment score: 5.77)
extracellular matrix organization
extracellular structure organization
collagen fibril organization
Interactions and functions of key modulator genes
Top 6 clusters of GO terms enriched in the 973 modulator genes
Cluster 1 (enrichment score: 5.59)
T cell activation
Cluster 2 (enrichment score: 5.29)
response to wounding
Cluster 3 (enrichment score: 3.60)
Cluster 4 (enrichment score: 3.39)
positive regulation of alpha-beta T cell activation
regulation of alpha-beta T cell activation
positive regulation of immune system process
regulation of lymphocyte activation
positive regulation of leukocyte proliferation
Cluster 5 (enrichment score: 3.28)
regulation of inflammatory response
positive regulation of response to stimulus
positive regulation of inflammatory response
regulation of response to external stimulus
positive regulation of response to external stimulus
Cluster 6 (enrichment score: 3.07)
response to organic substance
response to estrogen stimulus
response to steroid hormone stimulus
response to hormone stimulus
response to endogenous stimulus
Among the top GO clusters we also identified crucial tumor-related functions, such as cell adhesion (3rd cluster, score = 3.60, 66 modulators) and response to estrogen stimulus (6th cluster, score = 3.07, 68 modulators) (Table 3). The former is associated with metastasis and survival of breast cancer, while the latter is related to routine functions of hormonal receptors that were also seen in a previous co-modulation study . Interestingly, in the cluster of response to estrogen stimulus, in addition to ESR1 we identified another hormone receptor, androgen receptor (AR). Taken together with the well-studied role of ER as a modulator gene in breast cancer, our data showed that its functions, especially in the response to estrogen, are co-performed by other modulator genes, highlighting the essential involvement of modulation in such functions.
Limitations and future work
We measured the modularity of each putative modulator at two layers of interaction networks, one focusing on global changes and the other on a core subset of network. Four modularity parameters were designed accordingly, of which validity was confirmed by three independent datasets. However, for the nature of gene modulation as an indirect and complex mechanism, there may exist other parameters that could better measure modularity when cooperatively considered with the proposed four parameters. Furthermore, since the statistical features of the four parameters have not been characterized, we employed random permutations to assess the statistical significance, which limits the computation efficiency and statistical stringency. Out of simplicity, we compared the interaction networks formed in m-on and -off samples. However, modulated gene pairs of which correlation changes gradually with the continuous-state expression of m  may be omitted. Besides, in the study we assumed modulation effects to be independent events. Though, biological intuition is that several modulators may jointly modulate a common pair of genes , and pairs of genes modulated by a modulator may have competing effects against each other . Future investigation addressing the limitations may further unveil a comprehensive map of modulated gene interactions in cancers and other diseases.
This study addresses the need for a genome-wide screening for key modulator genes of gene interaction. We developed a systematic workflow that incorporates a correlation-based modulation analysis of gene interaction networks. About one thousand key modulators were identified, including the best-known modulator ESR1 and other novel ones, and validated in independent cohorts. These modulators were associated with hormone signaling and immune/defense-related and tumor-associated functions. Overall, this study is, to our knowledge, the first to screen for and investigate modulator genes in breast cancer on a genome-wide scale. The proposed workflow is widely applicable to other cancers and expected to unveil the landscape of modulated gene interactions.
The authors greatly appreciate the brilliant and constructive inputs from reviewers and participants of the International Conference on Intelligent Biology and Medicine (ICIBM 2016).
This research and this article’s publication costs were supported by the National Health Research Institutes of Taiwan (NHRI-EX106-10419BI). This research was also supported by the National Cancer Institute (1R01CA152063–02 and U54 CA113001–10) and Greehey Children’s Cancer Research Institute (GCCRI) intramural research fund. The funding sources had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
About this supplement
This article has been published as part of BMC Genomics Volume 18 Supplement 6, 2017: Selected articles from the International Conference on Intelligent Biology and Medicine (ICIBM) 2016: genomics. The full contents of the supplement are available online at https://bmcgenomics.biomedcentral.com/articles/supplements/volume-18-supplement-6.
YCC, LJW, THH, EYC, and YC conceived the study. YCC and LJW designed the study and performed data analysis. YCC, LJW, and YC interpreted the data. All authors wrote the manuscript. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Not applicable. All gene expression datasets used in this study were generated from previous studies and are publicly available at The Cancer Genome Atlas (TCGA; https://cancergenome.nih.gov/) and Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.
- Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Li KC. Genome-wide coexpression dynamics: theory and application. Proc Natl Acad Sci U S A. 2002;99(26):16875–80.View ArticlePubMedPubMed CentralGoogle Scholar
- Ideker T, Krogan NJ. Differential network biology. Mol Syst Biol. 2012;8:565.View ArticlePubMedPubMed CentralGoogle Scholar
- Flores M, Hsiao TH, Chiu YC, Chuang EY, Huang Y, Chen Y. Gene regulation, modulation, and their applications in gene expression data analysis. Adv Bioinforma. 2013;2013:360678.View ArticleGoogle Scholar
- Wilson CA, Dering J. Recent translational research: microarray expression profiling of breast cancer--beyond classification and prognostic markers? Breast Cancer Res. 2004;6(5):192–200.View ArticlePubMedPubMed CentralGoogle Scholar
- Shen C, Huang Y, Liu Y, Wang G, Zhao Y, Wang Z, Teng M, Wang Y, Flockhart DA, Skaar TC, et al. A modulated empirical Bayes model for identifying topological and temporal estrogen receptor alpha regulatory networks in breast cancer. BMC Syst Biol. 2011;5:67.View ArticlePubMedPubMed CentralGoogle Scholar
- Hsiao TH, Chiu YC, Hsu PY, Lu TP, Lai LC, Tsai MH, Huang TH, Chuang EY, Chen Y. Differential network analysis reveals the genome-wide landscape of estrogen receptor modulation in hormonal cancers. Sci Rep. 2016;6:23035.View ArticlePubMedPubMed CentralGoogle Scholar
- Chiu YC, Wu CT, Hsiao TH, Lai YP, Hsiao C, Chen Y, Chuang EY. Co-modulation analysis of gene regulation in breast cancer reveals complex interplay between ESR1 and ERBB2 genes. BMC Genomics. 2015;16(Suppl 7):S19.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang K, Saito M, Bisikirska BC, Alvarez MJ, Lim WK, Rajbhandari P, Shen Q, Nemenman I, Basso K, Margolin AA, et al. Genome-wide identification of post-translational modulators of transcription factor activity in human B cells. Nat Biotechnol. 2009;27(9):829–39.View ArticlePubMedPubMed CentralGoogle Scholar
- Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Holko M, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res. 2013;41(Database issue):D991–5.View ArticlePubMedGoogle Scholar
- Wang Y, Klijn JG, Zhang Y, Sieuwerts AM, Look MP, Yang F, Talantov D, Timmermans M, Meijer-van Gelder ME, Yu J, et al. Gene-expression profiles to predict distant metastasis of lymph-node-negative primary breast cancer. Lancet. 2005;365(9460):671–9.View ArticlePubMedGoogle Scholar
- Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, Nordgren H, Farmer P, Praz V, Haibe-Kains B, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. 2006;98(4):262–72.View ArticlePubMedGoogle Scholar
- Ivshina AV, George J, Senko O, Mow B, Putti TC, Smeds J, Lindahl T, Pawitan Y, Hall P, Nordgren H, et al. Genetic reclassification of histologic grade delineates new clinical subtypes of breast cancer. Cancer Res. 2006;66(21):10292–301.View ArticlePubMedGoogle Scholar
- Cancer Genome Atlas N. Comprehensive molecular portraits of human breast tumours. Nature. 2012;490(7418):61–70.View ArticleGoogle Scholar
- Ciriello G, Gatza ML, Beck AH, Wilkerson MD, Rhie SK, Pastore A, Zhang H, McLellan M, Yau C, Kandoth C, et al. Comprehensive Molecular Portraits of Invasive Lobular Breast Cancer. Cell. 2015;163(2):506–19.View ArticlePubMedPubMed CentralGoogle Scholar
- Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.View ArticlePubMedGoogle Scholar
- Huang da W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13.View ArticlePubMedGoogle Scholar
- Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.View ArticlePubMedPubMed CentralGoogle Scholar
- Kakinuma N, Zhu Y, Wang Y, Roy BC, Kiyama R. Kank proteins: structure, functions and diseases. Cell Mol Life Sci. 2009;66(16):2651–9.View ArticlePubMedGoogle Scholar
- Klopocki E, Kristiansen G, Wild PJ, Klaman I, Castanos-Velez E, Singer G, Stohr R, Simon R, Sauter G, Leibiger H, et al. Loss of SFRP1 is associated with breast cancer progression and poor prognosis in early stage tumors. Int J Oncol. 2004;25(3):641–9.PubMedGoogle Scholar
- Veeck J, Niederacher D, An H, Klopocki E, Wiesmann F, Betz B, Galm O, Camara O, Durst M, Kristiansen G, et al. Aberrant methylation of the Wnt antagonist SFRP1 in breast cancer is associated with unfavourable prognosis. Oncogene. 2006;25(24):3479–88.View ArticlePubMedGoogle Scholar
- Mohammed Ael S, Eguchi H, Wada S, Koyama N, Shimizu M, Otani K, Ohtaki M, Tanimoto K, Hiyama K, Gaber MS, et al. TMEM158 and FBLP1 as novel marker genes of cisplatin sensitivity in non-small cell lung cancer cells. Exp Lung Res. 2012;38(9–10):463–74.View ArticlePubMedGoogle Scholar
- Cheng Z, Guo J, Chen L, Luo N, Yang W, Qu X. Overexpression of TMEM158 contributes to ovarian carcinogenesis. J Exp Clin Cancer Res. 2015;34:75.View ArticlePubMedPubMed CentralGoogle Scholar
- Wu W, Tian Y, Wan H, Song Y, Li J, Zhang L. The expressions of Wnt/beta-catenin pathway-related components in brainstem gliomas. Can J Neurol Sci. 2013;40(3):355–60.View ArticlePubMedGoogle Scholar
- Kang P, Wan M, Huang P, Li C, Wang Z, Zhong X, Hu Z, Tai S, Cui Y. The Wnt antagonist sFRP1 as a favorable prognosticator in human biliary tract carcinoma. PLoS One. 2014;9(3):e90308.View ArticlePubMedPubMed CentralGoogle Scholar
- Zheng L, Sun D, Fan W, Zhang Z, Li Q, Jiang T. Diagnostic value of SFRP1 as a favorable predictive and prognostic biomarker in patients with prostate cancer. PLoS One. 2015;10(2):e0118276.View ArticlePubMedPubMed CentralGoogle Scholar
- Chang L, Lei X, Qin YU, Zeng G, Zhang X, Jin H, Wang C, Wang X, Su J. Expression and prognostic value of SFRP1 and beta-catenin in patients with glioblastoma. Oncol Lett. 2016;11(1):69–74.PubMedGoogle Scholar
- Ma XJ, Dahiya S, Richardson E, Erlander M, Sgroi DC. Gene expression profiling of the tumor microenvironment during breast cancer progression. Breast Cancer Res. 2009;11(1):R7.View ArticlePubMedPubMed CentralGoogle Scholar
- Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discovery. 2012;2(5):401–4.View ArticlePubMedGoogle Scholar
- Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, Jacobsen A, Sinha R, Larsson E, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6(269):pl1.View ArticlePubMedPubMed CentralGoogle Scholar
- Velasco-Miguel S, Buckbinder L, Jean P, Gelbert L, Talbott R, Laidlaw J, Seizinger B, Kley N. PA26, a novel target of the p53 tumor suppressor and member of the GADD family of DNA damage and growth arrest inducible genes. Oncogene. 1999;18(1):127–37.View ArticlePubMedGoogle Scholar
- Zhang WJ, Li BH, Yang XZ, Li PD, Yuan Q, Liu XH, Xu SB, Zhang Y, Yuan J, Gerhard GS, et al. IL-4-induced Stat6 activities affect apoptosis and gene expression in breast cancer cells. Cytokine. 2008;42(1):39–47.View ArticlePubMedGoogle Scholar
- Elhassan MO, Christie J, Duxbury MS. Homo sapiens systemic RNA interference-defective-1 transmembrane family member 1 (SIDT1) protein mediates contact-dependent small RNA transfer and microRNA-21-driven chemoresistance. J Biol Chem. 2012;287(8):5267–77.View ArticlePubMedGoogle Scholar
- Stanton SE, Adams S, Disis ML. Variation in the Incidence and Magnitude of Tumor-Infiltrating Lymphocytes in Breast Cancer Subtypes: A Systematic Review. JAMA Oncol. 2016;Google Scholar
- Savas P, Salgado R, Denkert C, Sotiriou C, Darcy PK, Smyth MJ, Loi S. Clinical relevance of host immunity in breast cancer: from TILs to the clinic. Nat Rev Clin Oncol. 2016;13(4):228–41.View ArticlePubMedGoogle Scholar
- Emens LA, Braiteh FS, Cassier P, DeLord J-P, Eder JP, Shen X, Xiao Y, Wang Y, Hegde PS, Chen DS, et al. Abstract PD1-6: Inhibition of PD-L1 by MPDL3280A leads to clinical activity in patients with metastatic triple-negative breast cancer. Cancer Res. 2015;75(9 Supplement):PD1-6-PD1-6.View ArticleGoogle Scholar
- Nanda R, Chow LQ, Dees EC, Berger R, Gupta S, Geva R, Pusztai L, Dolled-Filhart M, Emancipator K, Gonzalez EJ, et al. Abstract S1-09: A phase Ib study of pembrolizumab (MK-3475) in patients with advanced triple-negative breast cancer. Cancer Res. 2015;75(9 Supplement):S1-09-S01-09.View ArticleGoogle Scholar
- Chiu YC, Hsiao TH, Wang LJ, Chen Y, Chuang EY: Analyzing Differential Regulatory Networks Modulated by Continuous-State Genomic Features in Glioblastoma Multiforme. IEEE/ACM transactions on computational biology and bioinformatics/IEEE, ACM 2016.Google Scholar
- Chiu Y-C, Wang L-J, Lu T-P, Hsiao T-H, Chuang EY, Chen Y. Differential correlation analysis of glioblastoma reveals immune ceRNA interactions predictive of patient survival. BMC Bioinform. 2017;18(1):132.View ArticleGoogle Scholar