Transcriptional Networks Inferred from Molecular Signatures of Breast Cancer

Promoter Analysis

Hierarchical clustering data2 were used to construct an initial list of genes contained within the boundaries of the positive peaks of expression delimited by the gene order and clusters defined in Sorlie and colleagues2 (Figure 1, A and B). To do this the median expression of each of the 534 genes in each sample cluster was determined. This produced five median expression profiles for each of the 534 genes of each cluster (Figure 1B). The genes within the major contiguous positive peaks for each cluster were then broadly selected (color coded peaks in Figure 1B). This resulted in an initial overlapping list of genes containing: 136 genes from the basal subtype, 15 genes from the ERBB2⁺ subtype, 243 genes from the luminal A subtype, 217 genes from the luminal B subtype, and 108 genes from the normal-like subtype. These gene lists were then refined to remove overlapping and other noninformative genes by significance ranking for subtype discriminators using one-way analysis of variance and selecting genes with P values less than 0.01 (before correction for multiple comparisons) (see Supplementary Table S1 at http://ajp.amjpathol.org). The result was a final list of 221 genes: 95 unique genes for luminal A subtype, 21 unique genes for luminal B subtype, 66 unique genes for basal subtype, six unique genes for ERBB2⁺ tumor subtype, and 13 unique genes for the normal-like tumor subtype (see Supplementary Table S1 at http://ajp.amjpathol.org). Promoter sequences from three 600-bp promoter regions (proximal: −500 bp to +100 bp; upstream: −1100 bp to −500; and downstream: +100 bp to +700 bp, all relative to transcription start site) were retrieved for each gene from the five different tumor subtype-specific gene lists using the ProSpector web-based promoter annotation tool.10

The promoter regions were analyzed for matches to 409 position weight matrices using the MatInpector module of the GEMS LauncheR 4.1 (Genomatix, Munich, Germany). The analysis was performed with an optimized matrix threshold and a threshold of 0.75 for the core similarity. A P value for each gene in each gene list (cluster) was obtained by comparing the number of matches per 1-kb promoter region in each of the tumor subtypes against the number of matches per 1-kb promoter region in a reference background model using a complemented Poisson distribution (pdtrc) in the Perl math (math-cephes) library. The reference background model was obtained by extracting all unique RefSeq IDs (24,704) from the UCSC genome browser (build hg17) and mapping them to 15,318 IDs using Prospector followed by promoter TFBS annotation using MatInspector as described above. A Perl script was used to calculate P values for each of the 409 position weight matrices in each of the gene lists. The Perl script processed an input file containing a list of the 409 matrices and the number of matches for each of the matrices in a particular gene list and a file containing a list of the 409 matrices and the number of matches for each of the matrices in the reference background model. The script produced a file containing the list of the 409 matrices and the P values associated with each of the matrices calculated using the complemented Poisson distribution from the Perl math (math-cephes) library.

Statistical Analysis

One-way analysis of variance, hierarchical clustering, PCA, and intensity plots were generated using Partek Pro 5.1. Principal component (PC) loading correlation values were calculated with Partek Pro 5.1. PCA biplot analysis was performed as previously described.11 For this analysis, TFBS matrices with a correlation value less than 0.75 in any of the first four PCs were removed resulting in 208 matrices. This list was filtered further for a P value ≤0.05 in one of any of the five subtypes leaving 44 matrices for Biplot analysis. Randomization was performed by taking 40,000 random gene list selections of 6, 13, 21, 66, and 95 genes from the reference background list of 15,318 RefSeq genes and analyzing for frequency of the 409 matrices in each of the 40,000 random iterations for each of the five gene list sizes using MatInspector. Perl scripts were used to create the 40,000 random gene lists, to calculate matches for each of the 409 matrices in each of the random gene lists, and to calculate P values for each of the 409 matrices in each of the random gene lists.

Pathway and Network Analysis

Gene lists were analyzed with the Ingenuity Pathway Analysis software (Ingenuity Systems, Redwood City, CA). Networks were constructed by overlaying the genes in the gene list, called Focus Genes, onto a global molecular network developed from information contained in the Ingenuity Pathways knowledge base. Networks of these focus genes were then algorithmically generated based on their connectivity. A network is a graphical representation of the molecular relationships between genes. Genes are represented as nodes, and the biological relationship between two nodes is represented as an edge (line). All edges are supported by at least one reference from the literature, from a textbook, or from canonical information stored in the Ingenuity Pathways knowledge base. P values for the enrichment of canonical pathways were generated based on the hypergeometric distribution and calculated with the right-tailed Fisher’s exact t-test for 2 × 2 contingency tables. The composite gene list constructed using TFBSs included all known genes cognate for the particular binding site.

Chromatin Immunoprecipitation

MCF-7 and MDA-MB-231 breast cancer cell lines (gifts from Dr. Alfred Johnson and Dr. Ira Pastan, respectively, from the National Cancer Institute, Bethesda MD) were grown in Dulbecco’s modified Eagle’s medium and 10% fetal calf serum at 37°C and 5% CO_2. The cells were serum-starved overnight before a 1-hour stimulation with 12 ng/ml of EGF (Peprotech, Rocky Hill, NJ). Cells were then harvested, cross-linked with formalin, and chromatin immunoprecipitation was performed as previously described.12 The antibodies used were a 1:1 mixture containing 5 μg each of affinity-purified antibodies against p65 (sc-109; Santa Cruz Biotechnology, Santa Cruz, CA) and c-rel/Rel.13 Primers used for the ERBB2 promoter were 5′-TATTTTATCCTTGGTGTCGTGGCAGC-3′ and 5′-CATTGGCTGGCACTGGTCCC-3′.

Derivation of Breast Cancer Subtype-Specific Gene Signatures

Breast cancer subtype-specific signatures were derived from a previously published hierarchical clustering of 534 genes (represented by 552 clones; 500 of these correspond to a single unique UniGene cluster) expressed in 115 breast cancer samples (Figure 1A).2 As noted by Sorlie and colleagues,2 specific core genes within these 534 genes showed particular discriminatory power. For our promoter analysis we sought to expand this core list of subtype-specific genes. Therefore lists of up-regulated genes, comprising the core genes grouped with the adjacent members within each cluster, as delimited by the dendrograms for each subtype (Figure 1A), were constructed based on the peak median gene expression for the 534 genes (Figure 1B). Although focusing exclusively on up-regulated expression creates a bias toward finding presumed positive correlations between TFBS enrichment and gene expression, this simplified approach results in minimizing noise that would result from attempts to distinguish and interpret both positive and negative discriminators during the promoter analysis. Based on this assumption, a preliminary list of subtype-specific genes was constructed from the peak of median expression for each subtype. To reduce stochastic noise in these lists we performed a one-way analysis of variance analysis to remove noninformative genes (P value >0.01, before correction for multiple comparisons) that did not contribute significantly to the classification of the five tumor subtypes (see Supplementary Table S1 at http://ajp.amjpathol.org). This reduced the original 534 genes to a total 201 RefSeq mapped genes comprising five nonoverlapping subtype-specific genetic signatures (95 for luminal A, 21 for luminal B, 66 for basal, 6 for ERBB2⁺, and 13 for normal-like subtype; see Supplementary Table S1 at http://ajp.amjpathol.org). As shown in Figure 1C, hierarchical clustering of this refined list resembles the basic partitioning of the tumor subtypes shown in the original hierarchical clustering analysis.2,3 A similar approach was recently described by Muggerud and colleagues.14

Promoter Analysis Reveals Tumor Subtype-Specific Enrichment of TFBSs

To search for potential common regulatory pathways associated with the five subtype-specific genetic signatures, the promoter regions of each gene were extracted from −500 bp upstream to +100 bp downstream of the transcription start site (see Materials and Methods). These promoter regions were then scored for matches to 409 different TFBS matrices using the MatInspector module of GEMS Launcher 4.1 (Genomatix). The observed frequency of the 409 matrices in the promoters of each subtype-specific signature was then compared to the expected frequency using a reference background model of 15,318 human RefSeq genes. The null hypothesis that the observed frequency of TFBS in the selected upstream sequences could be explained as a random fluctuation, as compared to the background frequency, was used to estimate a P value for significant enrichment. A complemented Poisson distribution was used to model the random processes governing the frequency of TFBS in the human genome. We verified this theoretical assumption by comparing the P values derived from the Poisson distribution to frequencies simulated in a random permutation test with 40,000 randomly selected lists of genes of the same size as the original lists. This analysis showed good agreement between the analytic and experimental distributions for each of the gene signatures (see Materials and Methods and Supplementary Figure S1 at http://ajp.amjpathol.org). The P value for each of the 409 matrices in each of the tumor subgroup gene lists was then analyzed by hierarchical clustering after −log₂ transformation (Figure 2). Clustering of the matrices according to tumor subtype reveals multiple distinct groups of TFBSs that are enriched significantly and selectively within the five subtypes (Figure 2, Table 1). These patterns of TFBS enrichment are highly position-dependent because they are not recapitulated at contiguous 600-bp regions downstream or upstream of the proximal −500 to +100 regions characterized in this study (see Supplementary Figure S3 at http://ajp.amjpathol.org).

Many of the TFBSs enriched in the specific tumor subtypes (P value ≤0.05) bind factors that participate in well known transcriptional pathways. For example, binding sites for NF-κB (NF-κB) (P ≤ 0.042), E2F factors (P ≤ 0.01), EGR1 (early growth response protein 1) (P ≤ 0.01), and SMAD (P ≤ 0.01) are enriched in the ERBB2⁺ group indicating that signaling through NF-κB-, E2F-, EGR-, and transforming growth factor (TGF)-β-regulated pathways are likely to play functional roles in the biology of tumors classified by this molecular signature.15,16,17 The two groups of breast cancer that showed the greatest overlap in promoter composition were luminal B and basal-like subtype. A common feature of these groups is the relatively high GC content of many of the consensus sequences representing the matrices enriched in both groups despite the fact that the average GC content of the promoter regions of these groups are not substantially different from the other subtypes (Table 1 and Supplementary Table S1 at http://ajp.amjpathol.org). These include binding sites for early growth response factor-3, EGR3 (P = 0.0008/luminal B; P = 0.007/basal), and the nuclear respiratory factor NRF1 (P = 0.0001/luminal B; P value = 0.000006/basal).18 Distinguishing elements for luminal B include binding sites for the xenobiotic responsive aryl-hydrocarbon factor ARH/ARNT (P = 0.002) that plays a prominent role in estrogen metabolism and cross-talks significantly with ER in breast cancer cells,19 and a binding site for the Nurr1/Nur77 class orphan receptors NBRE (P = 0.009) which have differential roles in controlling cell death and proliferation.20 Elements that distinguish the basal-like group include multiple members of the PAX family of transcription factors known to play a role in cellular differentiation, stem cell biology, and cancer. These include PAX3 (P = 0.014), PAX9 (P = 0.0009), and PAX4/6 (PAX4_PD) (P = 0.047).21 Other highly significant and distinguishing pathways for the basal-like phenotype include the well known oncogene cMyc (P = 0.0008), a pathway commonly amplified in ER-negative breast cancer,22 and the BRN family of pou domain factors (P = 0.0001), which recently have been implicated to have a role in angiogenesis and cancer.23

Interestingly, the promoter signatures derived from this analysis show relationships among the different tumor subtypes that are different from those implied from the hierarchical clustering of the original gene expression data (Figure 1A). Dendrograms of the gene expression signature suggest a close relationship between the basal-like and ERBB2⁺ phenotypes whereas clustering of the regulatory signatures of these gene groups suggest a closer relationship between luminal B and the basal-like molecular phenotypes (Figure 2G). Although this may reflect differences in the clustering parameters (also note that the rederived gene list of 201 genes in Figure 1C has a slightly different relationship than that of Figure 1A), this finding suggests that such similarities in the molecular phenotypes may be derived from very different signaling pathways and may therefore implicate different regulatory and/or oncogenic origins despite the fact that they are both associated with a poor prognosis.

PCA of Tumor Subtype-Specific Regulatory Signatures Reveals Distinct Regulatory Trends and Minimizes Redundancies

In Figure 3A PCA was used to reduce the initial 409 variables or dimensions, represented by the 409 TFBS matrices, to three dimensions that capture most of the trends (variance) in the data set. As shown in Figure 3A, the tumor subgroups are separated in a three-dimensional space defined by the first three principal components (PCs) of the PCA representing 80.4% of the variance of the data set. Because greater than 80% of the trends in the data can be described in the three PCs, we focused on those patterns that were most apparent in three dimensions. Each PC is a composite derived from specific contribution from each of the original 409 matrices. Thus the relative position of the data points representing the five tumor subgroups are compared in three dimensional space based on their aggregate enrichment for the 409 TFBS. The distance between the data points for each tumor subtype in the three dimensions represents the relative similarity of their regulatory signatures. In other words, shorter distances between signatures indicate greater similarity in promoter composition. As assessed by this method, the tumor subtypes have regulatory signatures that are very distinct. Consistent with the hierarchical clustering in Figure 2, the greatest similarity exists between the basal-like and luminal B regulatory signatures.

Figure 3B shows that the remaining trends in the dataset are captured by the fourth PC. Because 100% of the trends in the data can be explained by four composite dimensions, the relationship between the tumor subtypes were compared using multiple two dimensional plots to extract any additional information that may be inferred from the remaining PC (Figure 3C). Here the basic relationship between the subgroups shown in Figure 3A persists except for PC1 versus PC4, which demonstrates a greater separation between the basal-like and luminal B signature along PC4. This suggests that matrices contributing strongly to PC4 are significant discriminators for the basal-like versus luminal B regulatory signatures. To identify which matrices contributed to this discrimination, we analyzed the principal component (PC) loading.11,24 PC loading provides the correlation of each of the original 409 variables (TFBS dimension) with the principal component vectors (see Supplementary Table S2 at http://ajp.amjpathol.org). Those matrices that make the strongest contribution to the composition of the PC will show significantly higher loading correlations L (eg, L ≥ 0.75). Major contributors to PC4 are BRN2 (L = −0.83), EVI1 (EVI1-myleoid transforming protein) (L = −0.90), NBRE (L = 0.80), NRSE (L = −0.82), PAX3 (L = −0.86), PAX9 (L = −0.79), SOX9 (L = −0.82), and ZNF35 (L = −0.83). All but NBRE have negative correlations, suggesting they are reciprocally enriched in one of the two subtypes. In this case all but NBRE are significantly enriched in the basal-like tumor subtype compared to luminal B. Interestingly most of the factors that bind these sites play substantial roles in embryogenesis and differentiation (see Discussion). Also both three- and two-dimensional PCAs show persistent separation of the luminal B, basal-like, and luminal A subtypes as a common group, distinct from the normal-like and ERBB2⁺ subtypes, along PC1. An examination of the PC loading of PC1 shows a strong correlation with AP2 (L = 0.987) suggesting that enrichment of this GC-rich matrix is a major distinction separating luminal B, basal-like, and luminal A from the normal-like and ERBB2⁺ regulatory signatures.

The parameters of PC loading can also be used to filter the original variables to produce a set of reduced size that best characterizes the trends in the regulatory signatures, highlighting those matrices that are significant discriminators despite borderline or low-ranking P values. To do this the 409 matrices were screened by PC loading for those matrices showing a PC loading ( L ) ≥0.75 in any one of the four PCs and a P value ≤0.05 in any one of the tumor subtypes. This reduced the list of 409 matrices to a final list of 44 (see Supplementary Table S3 at http://ajp.amjpathol.org). In Figure 3D these 44 original variables are superimposed on the PCA separation of the subtypes in the form of a biplot image. This projection provides a graphic illustration of how the 44 matrices discriminate the different tumor subtype signatures. For instance, AP1 is a dominating discriminating factor and pathway for the normal- like tumor subtype. The ERBB2⁺ signature projects in the direction of the NF-κB and E2F matrices. The basal-like subtype projects along various PAX matrices, whereas the luminal B subtype projects along the AHR/ARANT pathways. Finally luminal A appears to project in a direction bounded by STAF and NF1 pathways. As shown in Figure 3E, these 44 matrices produced a separation of breast cancer subtypes that is similar, although not identical, to the hierarchical clustering in Figure 2G.

Tumor Subtype-Specific Genes and Transcriptional Regulators Are Associated with Specific Regulatory Networks

To search for regulatory networks or pathways that appear most associated with the regulatory signatures found in the tumor subtype promoters, the genes of each respective subtype were combined with the transcription factor genes cognate for the most significantly enriched TFBSs within each subtype. This created a composite list of the regulated and regulator for each genetic signature.10 These enhanced lists were then used to interrogate the Ingenuity Knowledge database25 to construct regulatory networks based on curated interactions between the genes in each respective list.11 Figure 4 shows the five networks that were most populated by the five composite gene lists described above. The basal-like group is dominated by a highly connected Myc node. Luminal A is characterized by multiple hormone and orphan receptor nodes and, interestingly, luminal B contains a central hub with p53. The normal-like signature is dominated by AP-1 pathways and components. Most notably the ERBB2⁺ signature contains multiple nodes linking NF-κB, E2F, EGR, Forkhead, and SMAD pathways, and show many linkages to canonical pathways that regulate cell-cycle progression.

The canonical regulatory pathways associated with the composite genetic and regulatory signatures of each subgroup are shown in Table 2. As demonstrated in Figure 4 the ERBB2⁺ signature is dominated by cell-cycle regulation pathways (Fisher exact t-test P = 9.6, E-13) and NF-κB-associated signaling (P = 4.27, E-03). The basal-like subtype composite signature is also enriched in cell-cycle regulatory pathways in addition to WNT/catenin signaling and endoplasmic reticulum/ER-pathways. Normal-like breast cancer subtype shows relative overpopulation with diverse cytokine- and growth factor-responsive pathways downstream of Jak/Stat and MAP kinase signaling circuitry. The luminal A tumor subtype signature shows a high association with pathways linked to metabolism of hydrophobic amino acids, VEGF and IGF-1 signaling, although IGF-1 signaling is overrepresented in three of the five signatures. Finally, the luminal B composite signature is characterized by sonic hedgehog signaling, notch signaling, sterol synthesis, and p38 MAPK and cell-cycle signaling. The latter two are also shared by the basal-like subtype signature.

Inferred Identification of a Role for NF-κB and ERBB2 in a Self-Amplifying Regulatory Loop in Human Breast Cancer Cells

The selective enrichment of NF-κB TFBSs in the promoters of the ERBB2⁺ tumor subtype is consistent with observations that overexpression of ERBB2 is associated with increased activation of NF-κB.26,27,28,29,30 An intriguing aspect of this relationship is that in addition to ERRB2, the growth receptor adaptor protein, GRB7, a direct modulator of the ERBB2 receptor family31 (Figure 4D), is also a target of NF-κB pathways because both ERBB2 and GRB7 contain binding sites for NF-κB in their promoters (Table 1 and Supplementary Figure S4 at http://ajp.amjpathol.org). It is therefore very likely that these genes participate in a self-enhancing feed forward loop that amplifies NF-κB molecular signals driven by ERBB2 (see highlighted region in Figure 4D).

At high levels ERBB2 dimerizes with itself and becomes active in the absence of ligand.32 When present at physiological levels, ERBB2 readily heterodimerizes with all other members of the ERBB family, particularly EGFR, to produce ligand-specific complexes responsive to secreted EGF-1.33 To test whether or not an autoregulatory loop linking NF-κB to ERBB2 may exist in human breast cancer cells, we examined the in vivo association of NF-κB complexes with ERBB2 promoter before and after EGF-1 stimulation by chromatin immunoprecipitation (ChIP) in cells known to express normal (MCF-7) and amplified levels (MDA-MB-231) of ERBB2.33 An antibody cocktail containing affinity-purified antibodies specific for human p65/RelA and c-rel/Rel was used to perform ChIP in resting and EGF-1 stimulated MCF-7 and MDA-MB-231 breast cancer cells. As shown in Figure 5A, when normalized to nonspecific antibody and input DNA, there is a significant increase in NF-κB association with the ERBB2 promoter of MCF-7 cells after treatment with EGF-1. In contrast, NF-κB binding to ERRB2 in resting MDA-MB-231 is significantly higher than either resting or stimulated MCF-7. Moreover, the response to EGF-1 appears deregulated because the addition of EGF-1 fails to induce further NF-κB binding and instead shows some variable depression in MDA-MB-231. These novel data demonstrate that an autoregulatory loop, similar to what is schematically outlined in Figure 5B, exists in human breast cancer and is the first demonstration that NF-κB associates with the ERBB2 promoter in vivo in human breast-derived cells.