Integrative Analysis Reveals Subtype-Specific Regulatory Determinants in Triple Negative Breast Cancer

Abstract

Different breast cancer (BC) subtypes have unique gene expression patterns, but their regulatory mechanisms have yet to be fully elucidated. We hypothesized that the top upregulated (Yin) and downregulated (Yang) genes determine the fate of cancer cells. To reveal the regulatory determinants of these Yin and Yang genes in different BC subtypes, we developed a lasso regression model integrating DNA methylation (DM), copy number variation (CNV) and microRNA (miRNA) expression of 391 BC patients, coupled with miRNA-target interactions and transcription factor (TF) binding sites. A total of 25, 20, 15 and 24 key regulators were identified for luminal A, luminal B, Her2-enriched, and triple negative (TN) subtypes, respectively. Many of the 24 TN regulators were found to regulate the PPARA and FOXM1 pathways. The Yin Yang gene expression mean ratio (YMR) and combined risk score (CRS) signatures built with either the targets of or the TN regulators were associated with the BC patients' survival. Previously, we identified FOXM1 and PPARA as the top Yin and Yang pathways in TN, respectively. These two pathways and their regulators could be further explored experimentally, which might help to identify potential therapeutic targets for TN.

Keywords: FOXM1; PPARA; breast cancer; lasso; regulator; triple negative.

Figure 1

Study design. This study included three major parts. The first was to construct sample-specific models. Key regulators in different breast cancer subtypes were then identified from the models. Finally, we focused our follow-up exploration on the regulators in the TN subtype.

Figure 2

Top 30 over-represented gene ontology (GO) terms in biological process (BP) (GO-BP) terms in the Yin (a) and Yang (b) genes. “GeneRatio” in axis x is the ratio of the number of the identified Yin or Yang genes included in the corresponding term (“Count”) over the total number of the identified Yin or Yang genes. The y-axis is a description of GO-BP terms.

Figure 3

Model comparison. (a) The boxplot shows Spearman correlations between predicted and actual gene expression changes for all samples in models with different regulators. Using only miRNA-binding sites as features was significantly better than random; adding TF-binding sites significantly improved CV performance over using only miRNAs-binding sites; adding CNV significantly improved CV performance over using only TF- and miRNA-binding sites, while the model using TF/miRNA-binding sites, CNV and DM outperformed the others. (b) Boxplot shows Spearman correlations for all samples in Model 1 and Model 2. Model 1 was trained using CNV, DM, miRNA-binding sites and TF-binding sites from TRRUST, while Model 2 was trained using CNV, DM, miRNA-binding sites and TF-binding sites from ReMap. ^**** p-value < 2.2 × 10⁻¹⁶ by Wilcoxon signed rank test.

Figure 4

Error changes caused by regulators of the TN subtype. The x-axis indicates the miRNA/TF regulators, and the y-axis indicates increase in squared error across samples of the TN subtype after excluding the regulator from regression models. All regulators are ranked based on increase in squared error, and 24 key regulators (red dots) for the TN subtype were identified with FDR < 0.05.

Figure 5

TF-target and miRNA-target regulatory networks. (a) The network comprises 23 selected regulators (21 blue diamonds and 2 red diamonds) with nonzero predicted interactions with PPAR/FOXM1 pathway-related genes (circles) obtained from MSigDB. The red and purple circles indicate Yin (upregulated) and Yang (downregulated) genes, respectively. The color of edges ranging from grey to black corresponds to an order of increasing number of hits. (b) A subnetwork contains PPARA and PPARG, their targets, and regulators that regulate them. A 2-layer hierarchical structure forms with 10 upstream regulators, including 9 miRNAs and E2F1 on the top layer, and two downstream regulators PPARA and PPARG arranged at the bottom layer.

Figure 6

Survival analysis of the 31 gene YMR signature. TN patients were divided into high-risk (“High,” red line) and low-risk (“Low,” blue line) groups according their YMR signature scores. Survival fractions as a function of survival time (years) were then plotted for the two groups and the significant separation of the two curves were assessed by log-rank test. The YMR signature was tested using TNBC samples of TCGA and METABRIC datasets separately by the R package survcomp. The YMR signature can significantly stratify the 112 TCGA TNBC samples into high- and low-risk groups by both overall survival (OS) (a) and disease specific survival (DSS) (b). The YMR signature can also significantly stratify the 127 METABRIC TNBC samples into high- and low-risk groups by both OS (c) and DSS (d).

Figure 7

Kaplan–Meier survival analysis of the combined risk score. TN patients were divided into high-risk (“High,” red line) and low-risk (“Low,” blue line) groups according the 6 TF or 24 regulator combined risk scores. Survival fractions as a function of survival time (years) were then plotted for the two groups, and the significant separation of the two curves was assessed by log-rank test using the R package survcomp. The 6 TF combined risk score significantly stratified the TNBC samples from TCGA (a) and METABRIC (b) into high- and low-risk OS groups, respectively. The 24-regulator combined risk score significantly stratified the 112 TCGA TNBC samples into high- and low-risk groups in terms of OS (c) and DSS (d).