Identification of somatic mutation-driven enhancers and their clinical utility in breast cancer.

We generated an integrative pipeline to identify mutation-driven enhancers in breast cancer. First, we identified potential enhancers and their accessibility which correlated with their regulatory activity using assay for transposase-accessible chromatin with sequencing (ATAC-seq) data in breast cancer. Next, the breast cancer-associated gene mutation data obtained from The Cancer Genome Atlas (TCGA) revealed that these data involved 18,847 genes and 986 samples. Mutation enhancers were identified where at least one somatic mutation was located in enhancers. A total of 7635 mutant enhancers were obtained in breast cancer ( Figure 1 A). For example, mutant enhancers such as “chr10:8073699-8074200,” “chr10:8069545-8070046,” “chr12:57528013-57528514,” “chr13:44573793-44574294,” and “chr14:24162044-24162545” showed higher mutation frequencies in breast cancer. The mutation profiles of mutant enhancers were characterized, and it was found that missense mutations accounted for 40.8% of these enhancer mutations, and the number of missense mutations was significantly higher than other mutation classifications. The mutation types were mainly concentrated in SNPs (89.9%) with a small number of insertions (4.5%) and deletions (5.6%), most of which were C > T and C > G ( Figure 1 B). These mutation-driven enhancers are mainly located in exon, proximal promoter, and 3′-untranslated regions (UTRs) ( Figure S1 A). Enhancers differ between tumor subtypes, and mutations differ between subtypes. To better understand the link between mutations and subtypes, we categorized breast cancer samples into five subtypes: Basal, Her2, LumA, LumB, and normal. We analyzed the mutations in the different subtypes and found that the number of mutation-driven enhancers was highest in the LumB subtype, while the number of mutation-driven enhancers was lowest in the normal subtype. Notably, C > T mutations were the most common type of mutation in the LumB subtype ( Figure S1 B). Finally, mutant enhancers showing significant accessibility changes between samples with and without somatic mutations were identified as mutation-driven enhancers. We identified 135 mutation-driven enhancers that showed elevated chromatin accessibility in mutated breast cancer samples ( Figures 1 C and 1D; Table S2 ). We next characterized the genomic distribution of mutation-driven enhancers across different genomic features, including promoters, distal regions, coding exons UTRs, and intronic and exons regions. We found that 35.0% of mutation-driven enhancers mapped to promoters, 27.4% to exons, and 21.3% to introns ( Figure 1 E). For example, the C > T mutant enhancer (chr17: 39404925–39405426) was in the vicinity of the KRTAP9-4 promoter, and the mutant samples showed significantly higher chromatin accessibility compared to wild-type samples of the LumB subtype ( Figure 1 F). Figure 1 Identifying mutation-driven enhancers in breast cancer (A) Mutation spectrum of enhancer elements in breast cancer. (B) The distribution of mutation and mutation-driven enhancers in breast cancer. (C) The accessibility of mutation-driven enhancers which correlated with their regulatory activity using ATAC-seq data in breast cancer. (D) Mutation enhancer differential expression. (E) Percentage of distribution of mutation-driven enhancers within the genome. (F) The accessibility distribution of the mutation-driven enhancer (chr17:39404925-39405426) in the breast cancer. ATAC-seq data were visualized using Integrative Genomics Viewer tool (IGV). Identifying mutation-driven enhancers in breast cancer (A) Mutation spectrum of enhancer elements in breast cancer. (B) The distribution of mutation and mutation-driven enhancers in breast cancer. (C) The accessibility of mutation-driven enhancers which correlated with their regulatory activity using ATAC-seq data in breast cancer. (D) Mutation enhancer differential expression. (E) Percentage of distribution of mutation-driven enhancers within the genome. (F) The accessibility distribution of the mutation-driven enhancer (chr17:39404925-39405426) in the breast cancer. ATAC-seq data were visualized using Integrative Genomics Viewer tool (IGV). To characterize the molecular mechanism underlying how mutation-driven enhancers cause an increased risk of breast cancer, we identified potential target genes affected by mutation-driven enhancers based on guilt-by-association analysis 5 ( Table S1 ). We screened target genes/lncRNAs with significant effects of mutation-driven enhancers and constructed a mutant enhancer-gene regulatory network ( Figure 2 A). It contains 107 mutant enhancers and 201 downstream target genes (42 significantly upregulated, 159 significantly downregulated), and 47 lncRNAs (9 significantly upregulated, 38 significantly downregulated). To reveal the biological functions of the enhancer-gene regulatory network and gain a deeper understanding of the impact of mutant enhancers, gene ontology (GO) enrichment analysis were performed. The results of GO analysis were divided into biological processes (BP), cellular components (CC), and molecular functions (MF). The most enriched MF terms for downregulated target genes were protein heterodimerization activity. The most enriched BP terms for downregulated target genes were protein-DNA complex assembly and cellular hormone metabolic process. The most enriched CC terms for downregulated target genes were protein-DNA complex ( Figure 2 B). Enhancers could regulate gene expression through the recruitment of transcription factors (TFs) that bind to unique DNA recognition sequences. We therefore assessed whether the mutations located in enhancers damage the TF motifs. We employed motifbreakR 6 to measure the effects of somatic mutations on TF-binding motifs by using a highly efficient information content-based algorithm to discriminate between truly disruptive versus neutral variants. The stronger motifbreakR score indicates the bigger the difference between the wild-type and mutant allele, which assessed the variant effect on the analyzed TF motif ( Figure 2 C). Scores below 0.4 are marked as neutral, 0.7 as strong. 7 The results revealed that mutations in the enhancers strongly affect many TF motif families, many of which are known to play an important role in breast cancers. For example, the mutation in the “chr4:15066109-15066610” strongly affects the binding motifs of SOX family members (such as SOX2, SOX4, and SOX5) and the SoxF family (such as SOX7, SOX17, and SOX18). It had been reported that SOX2, SOX4, SOX5, and other SOX members are involved in tumorigenesis and the SoxF family could be attractive targets for anti-angiogenic therapy and the treatment of metastatic disease. 8 , 9 Mutations occurring in the “chr9:98056427-98056928” strongly affect the binding motifs of EGR family members (EGR1, EGR2, EGR3, and EGR4). Early growth response proteins (EGRs) have been reported to play a role in the prognosis of breast cancer. The higher the level of EGRs expression in breast cancer, the better are overall survival and disease-free survival. 10 , 11 Previous studies have confirmed that EGR2, EGR3, and EGR4 are key regulators of T cell activation in vitro and in vivo . 12 Mutations in the “chr1:81968091-81968592” strongly affect the androgen receptor (AR) binding motifs. Existing analyses have shown that AR plays a key role in the pathology and progression of breast cancer and is associated with all stages of breast cancer development. Mutations in the “chr17:38708828−38709329” strongly affect the binding motifs of PPARs family members, including PPARA, PPARB, and PPARG. PPAR signaling pathway is considered as actively involved in tumorigenesis. 13 Both PPARA and PPARG are associated with increased breast cancer risk. 14 There was a significantly differential expression of PPARA and PPARG in breast cancer (BRCA). The expression profiles of PPARA and PPARG was positively correlated with BRCA stromal scores. Figure 2 Mutation-driven enhancers contribute to dysregulation of potential target genes in breast cancer (A) The mutation-driven enhancer-target network. The fill colors of the nodes indicate the enhancers (yellow), PCGs (blue), mutations (red), and lncRNAs (green). The outer color of the node represents downregulated (gray) or upregulated (purple) genes. The genomic locations of enhancer elements are found in Table S1 . (B) Functional enrichment analysis of downstream target genes by GO. (C) The TF motif sequences of mutation-driven enhancers were analyzed by R package motifbreakR. Mutation-driven enhancers contribute to dysregulation of potential target genes in breast cancer (A) The mutation-driven enhancer-target network. The fill colors of the nodes indicate the enhancers (yellow), PCGs (blue), mutations (red), and lncRNAs (green). The outer color of the node represents downregulated (gray) or upregulated (purple) genes. The genomic locations of enhancer elements are found in Table S1 . (B) Functional enrichment analysis of downstream target genes by GO. (C) The TF motif sequences of mutation-driven enhancers were analyzed by R package motifbreakR. The 22% of potential target genes in the enhancer-gene regulatory network showed significant expression changes between breast cancer samples with and without the mutation-driven enhancer, including 42 upregulated protein-coding genes (PCGs), 159 downregulated PCGs, 9 upregulated lncRNAs, and 38 downregulated lncRNAs ( Figure 2 A). For example, the T > C mutation in the “chr8:109575180-109575681” had significantly upregulated the enhancer activity (log2FC = 2.90) in LumB subtype, which in turn was associated with the upregulation of potential target genes, such as ENY2, EBAG9, and SYBU (log2FC = 1.4) ( Figure 3 A; Table S1 ). The previous study showed that the estrogen receptor-binding fragment-associated antigen 9 (EBAG9) was identified as a potential prognostic factor for poor outcome of breast cancer patients with adjuvant tamoxifen therapy. EBAG9 immunoreactivity is a potential prognostic factor for poor outcome of breast cancer patients with adjuvant tamoxifen therapy. 15 It has been reported that SYBU was overexpressed in pancreatic adenocarcinoma, breast cancer, and bladder cancer. Recent research revealed that SYBU showed significantly different mutation frequencies in triple-negative breast cancer with high CD8 T cells and CD4 memory-activated T cells. 16 The DEL mutation raised enhancer activity (chr3:57896344-57896845; log2FC = 2.31) in LumB subtype, which in turn was associated with the downregulation of target genes ARF4, DENND6A, PXK, FLNB-AS1, FLNB, RPP14, SLMAP, PDE12, and ABHD6 ( Figure 3 B; Table S1 ). PXK was significantly associated with the prognosis of diffuse large B cell lymphoma. 17 The G > T mutation in the “chr9:65675405-65675906” significantly increased the activity of the enhancer (log2FC = 2.51) in LumB subtype, which in turn was associated with the downregulation of potential target genes, such as CBWD5 and ANKRD20A3P ( Figure 3 C; Table S1 ). The C > T mutation in “chr22:36194883-36195384” had the significantly upregulated the enhancer activity in LumA subtype, which in turn was associated with the upregulation of target genes TXN2, APOL1, APOL2, APOL3, APOL4, and EIF3D ( Figure 3 D; Table S1 ). Studies have shown that APOL1 is significantly associated with early recurrence in patients with BRCA. Increased expression of EIF3D can enhance the translation of U2SURP and play an important role in promoting triple-negative breast cancer (TNBC) cell tumorigenesis and metastasis in vivo and in vitro . 18 The mutation type C > T occurred on the “chr3:101859127-101859628” in LumB subtype, which increased the activity of the enhancer but downregulated the expression of its downstream target genes NFKBIZ, ZPLD1, ZBTB11, SENP7, CEP97, and NXPE3 ( Figure 3 E; Table S1 ). Target gene NFKBIZ is an independent prognostic feature of breast cancer. 19 Enhancer (chr12:11752019-11752520) also had the C > T mutation in Her2 subtype, and the enhancer activity was significantly upregulated (log2FC = 2.48); the expressions of downstream lncRNA AC007450.1 and AC078950.1 were significantly decreased. The long-stranded non-coding RNA AC078950.1 was thought to be significantly associated with the prognosis of patients with clear cell renal cell carcinoma. 20 The C > A mutation in the “chr16:31488177-31488678” significantly increased the activity of the enhancer (log2FC = 1.66) in LumB subtype, which in turn was associated with the upregulation of potential target genes, such as LINC02190, KAT8, and TRIM72. It had been reported that KAT8 was associated with a significantly high grade and poor prognosis of breast cancers. 21 Overexpression of TRIM72 is reported to inhibit tumor proliferation and invasion in vitro and in a xenograft tumor model in breast cancer progression. 22 The “DEL” mutation in the “chr3:57896344-57896845” significantly increased the activity of the enhancer (log2FC = 2.31) in LumB subtype, which in turn was associated with the downregulation of potential target genes, such as lncRNA FLNB-AS1, FLNB, and SLMAP. The lncRNA FLNB-AS1 was significantly downregulated (log2FC = −2.34) in breast cancer samples with the mutation-driven enhancer, and literature has confirmed that FLNB-AS1 is positively correlated with the survival probability of patients with breast cancer. 23 FLNB, which plays a causal role in the regulation of epithelial-mesenchymal transition (EMT), was found to inhibit breast cancer cell migration. 24 The T > A mutation in the enhancer (chr8:123182900-123183401) significantly decreased the activity of the enhancer (log2FC = 4.17) in Her2 subtype, which in turn was associated with the downregulation of potential target genes, such as lncRNA AC016405.3 , DERL1, and FAM83A. The expression of AC016405.3 was significantly downregulated (log2FC = −1.74) in breast cancer samples with the mutation-driven enhancer. AC016405.3 was confirmed to act as an oncogenic lncRNA and positively correlated with poor prognosis in breast cancer patients. 25 The expression of DERL1 has been reported to exhibit high accuracy in differentiating BRCA tissues from normal samples and its inhibition can partially abolish the oncogenic effects of TMEM63A on TNBC progression. 26 , 27 The previous study showed that FAM83A is a potential biomarker for breast cancer initiation and could promote metabolic activation in primary breast epithelial cells and cell proliferation in both primary and immortalized cells. 28 Altogether, mutation-driven enhancers might play important roles during breast cancer progression by dysregulation of their potential target genes. These results suggest that mutation-driven enhancers can be considered as potential risk factors in the breast cancer by disrupting a few important transcription factors and downstream target genes. Figure 3 Effects of mutant driver enhancers on potential target genes in different subtype (A‒E) Accessibility distribution of mutation-driven enhancers in different subtypes of breast cancer. ATAC-seq data were visualized (left panels) using Integrative Genomics Viewer tool (IGV). The distribution of chromatin accessibility of enhancers in mutant (red) and wild-type samples (blue; middle panels). The expression distribution of enhancer targets in mutant (red) and wild-type samples (blue; right panels). Effects of mutant driver enhancers on potential target genes in different subtype (A‒E) Accessibility distribution of mutation-driven enhancers in different subtypes of breast cancer. ATAC-seq data were visualized (left panels) using Integrative Genomics Viewer tool (IGV). The distribution of chromatin accessibility of enhancers in mutant (red) and wild-type samples (blue; middle panels). The expression distribution of enhancer targets in mutant (red) and wild-type samples (blue; right panels). To characterize the prognostic value of mutation-driven enhancers in breast cancer, we analyzed 135 mutation-driven enhancers expression in 1094 patients. Univariate Cox proportional hazards analysis identified the mutant enhancers that had an impact on the overall survival of BRCA patients, and p < 0.05 was set as the threshold of significance level, and two mutation-driven enhancers were significantly associated with overall survival. We then performed multivariate Cox proportional hazards model analysis on the expression of these 4 mutation-driven enhancers in relation to clinical parameters such as age at diagnosis, race, and pathological stage. We identified the 4 mutation-driven enhancers, including “chr17:7858334-7858835” (HR = 1.4, p = 0.017), “chr9:3346483-3346984” (HR = 2.6, p = 0.037), “chr19:12719195-12719696” (HR = 1.7, p = 0.0023), and “chr19:12834462-12834963” (HR = 1.5, p = 0.0033) as independent risk factors for BRCA prognosis. The Kaplan-Meier survival curve showed that mutations of the enhancers were associated with poorer prognosis (chr9:3346483-3346984, p = 0.041; chr17:7858334-7858835, p = 0.016; chr19:12719195-12719696, p = 0.0038; chr19:12834462-12834963, p = 0.0018; Figure 4 A). For example, in the Her2 subtype, G > C mutation of the enhancer (chr19:12719195-12719696) significantly downregulated the enhancer activity ( Figures 4 B and 4C), which in turn was associated with the dysregulation of potential target genes FBXW9, TRIR, and WDR83 ( Figures 4 B; Table S3 ). Previous studies have shown that WDR83 plays an important role in ERK signaling and tumorigenesis. 29 In the basal subtype, INS mutation significantly increased the chromatin accessibility of the enhancer (chr9:3346483-3346984; log2FC = 1.69, Figures 4 B and 4C), which may regulate the expression of the potential target gene AL354977.2 ( Table S3 ). In the LumB subtype, the chromatin accessibility of enhancers (chr19:12834462-12834963) was lower in mutant samples than in other samples ( Figures 4 B and 4C). The mutant enhancer was associated with downregulation of potential target genes AC010422.1 , KLF1, and RAD23A ( Table S3 ). It has been shown that KLF1 is strongly associated with anemia. 30 In the normal subtype, the G > A mutation significantly increased the chromatin accessibility of the enhancer (chr17:7858334-7858835, Figures 4 B and 4C), influencing its target genes ALOX15B, LINC00324, and MPDU1 ( Table S3 ). The ALOX15B (log2FC = −3.01) and LINC00324 (log2FC = −1.43) were significantly differentially expressed in the BRCA. The previous study showed that ALOX15B, LINC00324, and MPDU1 were linked to immune-mediated disease and cancers. 31 , 32 Figure 4 Characterization of the prognostic value of mutation-driven enhancers in breast cancer (A) Kaplan-Meier survival analysis of the mutation-driven enhancers. (B) ATAC-seq data of prognosis-related mutation-driven enhancers were visualized. (C) The distribution of chromatin accessibility of enhancers in mutant (red) and wild-type samples (blue; right panels). (D and E) Risk score of OS based on risk scores of the prognostic signature and ROC curves. Characterization of the prognostic value of mutation-driven enhancers in breast cancer (A) Kaplan-Meier survival analysis of the mutation-driven enhancers. (B) ATAC-seq data of prognosis-related mutation-driven enhancers were visualized. (C) The distribution of chromatin accessibility of enhancers in mutant (red) and wild-type samples (blue; right panels). (D and E) Risk score of OS based on risk scores of the prognostic signature and ROC curves. Differential transcription factor binding analysis showed that INS mutation in the enhancer (chr9:3346483-3346984) could affect the binding of transcription factors AR, NR3C1, NANOG, and progesterone receptor (PGR). AR is an important prognostic marker and potential therapeutic target in breast cancer. NANOG is an important stem cell transcription factor that is involved in the process of human carcinogenesis. 33 Increasing studies confirm that the PGR is significantly upregulated in breast cancer, which can act as a putative diagnostic biomarker in breast cancer detection. It has been reported that NR3C1 as an oncogenic factor participated in cell growth, glucocorticoid-induced apoptosis, inflammation, and differentiation. 34 NR3C1 expression was overexpressed in breast cancer, and the repression of NR3C1 blocked breast cancer cell migration and invasion. 35 The G > A mutation in the enhancer (chr17:7858334-7858835) could affect the binding of 15 transcription factors, such as POU6F2, POU2F1, EGR2, and EN1. Recent studies have found that POU6F2 plays an important role in the occurrence and development of pituitary adenomas and is a new target for the treatment of pituitary adenomas. 36 Studies have shown that EN1 transcription factor is a prognostic marker and potential therapeutic target for triple-negative breast cancer. 37 Some studies have suggested that early growth response 2 (EGR2) expression is correlated with apoptosis promotion, with the EMT process inhibitor. 38 EGR2 was downregulated in triple-negative breast cancer samples and involved in regulation for inducing EMT in highly invasive cells. The mutation-driven enhancer (chr17:7858334-7858835) was associated with the dysregulation of potential target genes, such as genes TNFSF12, TMEM102, CTC1, SAT2, CYB5D1, POLR2A, and lncRNA LINC00324. We then established the risk score consisting of 4 prognosis-related enhancers as a prognostic signature in BRCA. The risk score was calculated as follows: the activity of enhancer (chr17:7858334-7858835) ∗1.4 + the activity of enhancer (chr9:3346483-3346984) ∗2.6 + the activity of enhancer (chr19:12834462-12834963) ∗1.5 + the activity of enhancer (chr19:12719195-12719696) ∗1.7. The patients were divided into high-risk group and low-risk group according to the median of the risk score. A higher risk score resulted in shorter survival time and the death rate was higher in the high-risk group than in the low-risk group ( Figure 4 D). A time-dependent receiver operating characteristic (ROC) analysis was performed to calculate area under the curve (AUC) values to assess the predictive power of the model. The ROC curve indicated that the AUC values of the model were 0.835 ( Figure 4 E), indicating that the prognostic signaling properties of these two enhancers in breast cancer are highly reliable. When used alone as the diagnostic test, AUC for the 9 observers (including age, race, tumor stage, prior treatments, and four mutation-driven enhancers) ranged from 0.491 to 0.781, significantly less than 4-enhancer model. These results indicated that the four-mutation-driven enhancer signature demonstrated the better prognostic power in evaluating the prognosis of patients with BRCA. To characterize potential clinical applications of mutation-driven enhancers, we obtained the clinical data including clinical stage (Stage I, II, III, and IV), subtypes (Basal, Her2, LumA, Normal, and LumB), pathologic T stage, pathologic N stage, and pathologic M stage of breast cancer patients from TCGA. We characterized the differences in clinical features between mutated vs. non-mutated enhancers in the breast cancer samples. We found that prognostic mutation-driven enhancers were significantly associated with clinical features in breast cancer patients ( Figure S2 ). For example, stage II, T2 pathologic T stage, and Her2 subtype were highly frequent in breast cancer patients with mutated enhancer (chr17:7858334-7858835). Stage IV, N2 pathologic N stage, and M1 pathologic M stage were highly frequent in breast cancer patients with mutated enhancer (chr19:12834462-12834963). These results are consistent with poorer survival in breast cancer patients with the mutation of these enhancers. To characterize the clinically applicable therapeutic implications of mutation-driven enhancers, we identified 16 target genes of the enhancer (chr17: 7858334–7858835), one target gene of the enhancer (chr9:3346483-3346984), 14 target genes of the enhancer (chr19:12719195-12719696), and 14 target genes of enhancer (chr19:12834462-12834963). We performed the Gene Set Enrichment Analysis to calculate whether the target genes of mutation-driven enhancers were affected by drug perturbations based on the Broad Institute’s Connectivity Map. In total, we identified eight candidate small molecular drugs (PHA-00816795, zomepirac, colforsin, pinacidil, STOCK1N-35874, 0316684-0000, erastin, and PF-00539758-00) targeting the mutated enhancer (chr19:12719195-12719696) and two candidate drugs (aminoglutethimide and quinpirole) targeting the mutated enhancer (chr17: 7858334–7858835; FDR<0.01). Some potential compounds have been reported to demonstrate potential as an anticancer property. For example, erastin could inhibit the viability of breast cancer cells and induced breast cancer cell death. 39 , 40 Colforsin can inhibit tumor growth and tumor burden in an intraperitoneal model of high-grade serous ovarian cancer. 41 STOCK1N-35874 showed marked anticancer activity in preliminary clinical trials, and its analogs have been used in cancer chemotherapy. 42 Aminoglutethimide was the first inhibitor evaluated in clinical studies for treatment of hormone-dependent breast cancer. 43 These small molecular drugs could be potential candidates for enhancer-based cancer therapy.

Through integrated analysis of ATAC-seq, mutation data, gene expression data, and clinical data, we identified 7635 mutant enhancers in the breast cancer. The mutation types of mutant enhancers were mainly concentrated in SNPs (89.9%) with a small number of insertions (4.5%) and deletions (5.6%), most of which were C > T and C > G. Furthermore, we identified 135 mutation-driven enhancers that showed elevated chromatin accessibility in mutated breast cancer samples. We constructed an enhancer-gene regulatory network which consists of 248 regulatory relationships between 107 mutation-driven enhancers and their potential target genes, involving 201 PCGs and 47 lncRNAs. For example, the T > C mutation in the enhancer (chr8:109575180-109575681) had significantly upregulated the enhancer activity, which in turn was associated with the upregulation of potential target genes, such as ENY2, EBAG9, and SYBU. TF motif perturbation analysis revealed that somatic mutations in the enhancer strongly affect many TF motif families, many of which are known to play an important role in breast cancers. For example, the mutation in the enhancer (chr4:15066109-15066610) strongly affect the binding motifs of SOX family members (such as SOX2, SOX4, and SOX5). It had been reported that SOX family members are involved in tumorigenesis. To further the prognostic value of mutation-driven enhancers in breast cancer, univariate and multivariate Cox proportional hazards analysis were performed. We identified the 4 mutation-driven enhancers (“chr17:7858334-7858835,” “chr9:3346483-3346984,” “chr19:12719195-12719696,” and “chr19:12834462-12834963”) as independent risk factors for BRCA subtype prognosis. The G > A mutation significantly increased the activity of the enhancer (chr17:7858334-7858835) in the normal subtype. The INS mutation significantly increased the activity of the enhancer (chr9:3346483-3346984) in the basal subtype. The G > C mutation significantly increased the activity of the enhancer (chr19:12719195-12719696) in the Her2 subtype. The T>C mutation significantly increased the activity of the enhancer (chr19:12834462-12834963) in the LumB subtype. The higher activity of the four-mutation-driven enhancer signature was associated with poorer prognosis of BRCA. We identified the enhancer mutation rs763236375 as potential risk factors in breast cancer by influencing the expression of target gene GATA3. Gao et al. found that the SNP mutation in CACA at the rs763236375 site was the important reason for affecting GATA3 expression and high expression of GATA3 is beneficial to the prognosis of BRCA patients. 44 Similarly, we identified the enhancer mutation rs11554273 as a potential risk factor that may influence the progression of breast cancer. Zhao et al. demonstrated the rs11554273 as risk factors for fibrous dysplasia through mouse experiments. 45 The important role of the enhancer target genes identified in this study in breast cancer has been confirmed by many breast cancer-related studies. For example, the luciferase reporter assay confirmed that the overexpression of B cell translocation gene 2, a target gene of the enhancer (chr1:203305170-203305671), might inhibit proliferation, invasion, and induce apoptosis in breast cancer cells. 46 , 47 The results of in vitro and in vivo studies showed that knockdown of MEX3A, a target gene of enhancers (chr1:156241956–156242457 and chr1:156051355-156051856), inhibited the metastasis and impaired the stemness of breast cancer cells. 28 , 47 The breast tissue microarray analysis confirmed that enhancer (chr8:123182900-123183401) target gene FAM83A is a potential biomarker for breast cancer initiation. 28 , 48 RasGRP3 is the target gene of the enhancer (chr2:33252614-33253115). RasGRP3 knockdown in breast cancer cells decreased cell proliferation, induced apoptosis in MCF7 cells, and sensitized T-47D cells to the action of drugs tamoxifen and trastuzumab. 49 , 50 GATA3-AS1, a target gene of the mutation-driven enhancers (chr10:8073699–8074200 and chr10:8069545-8070046), contributed to TNBC progression and immune evasion through stabilizing PD-L1 protein and degrading GATA3 protein, offering a new target for the treatment of TNBC. 48 , 51 TONSL, a target of enhancer (chr8:144778224-144778725), is characterized by amplification, promotion of cell immortalization, and involvement in breast cancer tumorigenesis. 50 , 52 The high expression of BOP1, a target of enhancer (chr8:144414595-144415096) is associated with chemoresistance, poor prognosis, and tumor recurrence within 5 years in TNBC. 53 The expressions of TRAF4, a target gene of the enhancer (chr17:28741265-28741766), were enhanced in breast primary tumor using ELISA assay and immunohistochemical analysis. 52 When targeted with shRNA to knock down endogenous Rab34 which is target gene of the enhancer (chr17:28741265-28741766), it can significantly inhibit the migration and invasion ability of highly aggressive breast cancer cells. 51 , 54 Therefore, we identified many important mutation-driven enhancers and their target genes that may serve as potential risk factors for breast cancer progression. In summary, we present a systematical characterization of the mutation-driven enhancers and their clinical relevance in breast cancer. Our results underlined a significant association between the outcome of breast cancer patients and the mutation-driven enhancers, suggesting that mutation-driven enhancers could act as prognostic markers and potential therapeutic targets in human breast cancer. We used whole-exome sequencing data to characterize the functional effects of enhancer mutations, which mainly focus on gene proximal enhancers, such as exonic enhancers, intronic enhancer, and 5′- and 3′-UTR enhancers. However, the detection of mutations in intergenic enhancers is limited. As more large-scale whole-genome sequencing data of BRCA become available, it could further improve predictive capacities of our approach. Finally, the identification of subtype-specific prognostic genes and their functional role in BRCA subtypes need further investigation.

Breast cancer is the most common cancer among women today; the incidence of breast cancer is still increasing, and it remains a significant cause of premature death among women. With the development of next-generation sequencing technology, there is increasing evidence that somatic mutations are associated with an increased risk of breast cancer. Mutations in coding regions of the genome may alter nonfunctional proteins that cause cancer or other diseases. Noncoding mutations may act by disrupting (gaining or losing) transcription factor-binding motifs in promoters, enhancers, or silencers, resulting in altered target gene expression and increased cancer risk. 1 In breast cancer patients, especially triple-negative breast cancer, exon 4 and exon 3 of the TP53 gene are frequently mutated, and changes in this gene lead to the expression of various genes directly or indirectly under the transcriptional control of p53. These changes lead to malfunction of DNA damage repair pathways, cell-cycle arrest, chromatin remodeling, and apoptosis. 2 Germline PALB2 mutations have been reported to play an important role in hereditary breast cancers, with carriers of the mutation having a 5 times or greater risk of developing breast cancer than those with hereditary breast cancer. Mutations in enhancer elements may contribute to the development of diseases, with implications for downstream protein-coding genes. A study on leukemia found that genetic mutations in enhancer elements of PAX5 and TAL1 genes resulted in altered enhancer activity, which in turn drove abnormal expression of downstream target genes and long noncoding RNAs (lncRNAs), contributing to cancer development. Ronald L Chandler et al. found that the ARID1A mutation may inhibit the mechanism of invasive phenotype by antagonizing the P300 activity of the super enhancer, increasing the risk of endometriosis and malignant transformation. 3 It is well known that non-coding genomes contain a number of cis -regulatory elements, which are defined as binding sites of transcription factors and participate in transcription regulation by acting as promoters and enhancers of chromatin interactions. Somatic mutations affect transcription factor-binding motifs by affecting enhancer activity, thus interfering with gene expression. PAX6 is a major regulator of eye development. The mutation of PAX6 gene will affect the interaction with its upstream enhancer and lead to a wide range of clinical phenotypes, further affecting the binding with upstream transcription factors Sox2 and Sox3 and affecting tissue-specific transcription. 4 In this study, we systematically analyzed the somatic mutations of enhancers and its association with prognosis of breast cancer by integrating multi-omics data. We identified 7,653 mutant enhancers in breast cancer and 135 mutation-driven enhancers which showed significant chromatin accessibility changes between samples with and without somatic mutations. We constructed an enhancer-gene regulatory network which consists of 248 regulatory relationships between 107 mutation-driven enhancers and their potential target genes, involving 201 protein-coding genes and 47 lncRNAs. Transcription factors motif perturbation analysis revealed that mutations in the enhancer strongly affect many TF motif families, many of which are known to play an important role in breast cancers, such as SOX, PPARs, and EGR family. Finally, four mutation-driven enhancers were identified as independent risk factors for breast cancer subtype prognosis. Our results emphasized that somatic mutations drive the dysregulation of enhancer activity, thereby promoting breast cancer initiation and prognosis.

REAGENT or RESOURCE SOURCE IDENTIFIER Deposited data TCGA-BRCA mutations data NIH GDC data portal https://portal.gdc.cancer.gov/ BRCA enhancer data Juri Reimand et al. 55 https://doi.org/10.1016/j.molcel.2019.12.027 TCGA-BRCA gene expression data and clinical information NIH GDC data portal https://portal.gdc.cancer.gov/ Software and algorithms R (version 4.2.0) The R Foundation https://www.r-project.org/ maftools R package N/A ggplot2 R package N/A motifbreakR R package N/A clusterProfiler R package N/A survival R package N/A survminer R package N/A ggrisk R package N/A pROC R package N/A IGV (version 2.11.1) The IGV Foundation https://www.igv.org/ Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Li Wang ( [email protected] ). The study did not generate any new materials. • This study analyzes existing publicly available data. Relevant data information is in the key resources table . • This paper does not report original code. • Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. This study analyzes existing publicly available data. Relevant data information is in the key resources table . This paper does not report original code. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. The data analyzed in this study were obtained from The Cancer Genome Atlas (TCGA) database. All breast cancer candidate risk factor mutation enhancers and associated target genes in this study are summarized in Table S1 . Gene expression profiles and clinical survival information for BRCA were obtained from the TCGA database. The ATAC-seq raw data of 74 breast cancer patients were downloaded from the TCGA database, of which 66 breast cancer samples were tested twice. Data preprocessing and data comparison are performed based on the PEPATAC pipeline. FastQC performed quality control on the data, using the Bowtie algorithm to align ATAC-seq data to the human reference genome (hg19), and filtered reads to remove duplicates, unmapped or poor quality (Q < 30) Reads. Open chromatin regions were further identified using the MACS2 algorithm. The MACS2 algorithm parameters are as follows: --keep-dup all -B -noommodel--SPMR -q 0.05--slocal6250--llocal 6250. 56 , 57 Finally, 215,920 ATAC-seq enriched regions were identified in the breast cancer as potential breast cancer enhancers. Furthermore, to filter out peaks with low read counts, we converted the raw counts into CPM (counts per million) values to correct for the differences in library size between samples by using edgeR, which was used to characterize chromatin accessibility. In addition, mutation data, gene expression profiling data, and related clinical data of breast cancer patients were downloaded from the TCGA database. Somatic mutation data for breast cancer identified by the Mutect algorithm were obtained. Somatic mutations with low-read/low-bit gene fragments were identified through a Bayesian classifier while maintaining high specificity, assessing cross contamination between samples, and filtering the mutation set through FilterMutectCalls. Finally, the obtained high-quality mutation data are used in cancer-related research. The probes from TCGA whole exome sequencing (WES) data could contain non-coding regions, which are close to exon boundary such as 5′-untranslated region (UTR) and 3′-UTR and introns. Therefore, TCGA WES data was used to identify the mutations that affect the enhancer activity of gene proximal enhancers, such as exonic enhancers, intronic enhancer, 5′- and 3′-UTR enhancers. The mutation data involved 986 breast cancer samples containing 18,847 genes. The mammary gland mutation data was obtained, analyzed, and visualized using maftools 58 in the R package. After removing the expression data of low-quality samples and filtering out genes whose expression value is 0 in all samples, the expression profile data of 60,483 genes corresponding to 1,109 breast cancer samples and 113 normal samples were finally obtained. A permutation test was used to assess whether mutations in the enhancers significantly change their chromatin accessibility between samples with and without somatic mutations. Specifically, for a mutant enhancer, we compared fold change of chromatin accessibility in mutants and wild-type samples. To evaluate statistical significance, we randomly shuffled the sample labels (mutants and wild-type) and then re-calculated fold change of chromatin accessibility in mutant versus wild-type samples to build null distribution. After repeating the procedure 10,000 times, we would obtain a P-value of chromatin accessibility changes by dividing the times when the random fold change was greater than the observed of chromatin accessibility changes by 10,000. The mutant enhancers with P-value <0.05 were identified as mutation-driven enhancers in the breast cancer. Additionally, nonparametric tests were used to assess whether mutations in enhancers significantly altered their chromatin accessibility. Specifically, for a mutant enhancer, we compared the fold change in chromatin accessibility in mutant and wild-type samples. The rank sum (or rank mean) of each enhancer was calculated, and then the Wilcoxon rank-sum test was utilized to compare whether the difference between the mutant and non-mutant rank sums was significant. Mutant enhancers with a P-value <0.05 were identified as mutation-driven enhancers in breast cancer. We ultimately merged these two methods to obtain 135 enhancers with significant mutations. All possible interactions between the mutation-driven enhancer and its upstream and downstream target genes within 500kbp 5 were calculated by guilt-by-association analysis. Pearson correlation and significant p value between ATAC-seq peak activity and gene expression were calculated. Then using Chromosomal location information of enhancers and downstream target genes, using the 'intersect' function of the bedtools tool in the linux environment, parameters -c, -loj and -wo. By this method identified the interactions between enhancers in the upstream and downstream 500kbp, found significantly (threshold p < 0.05) associated co-expression regulatory pairs. We used motifbreakR 6 (v.2.6.1) to predict the effect of mutations in enhancer on transcription factor motifs. Positional weight matrices (PWMs) for transcription factors (TFs) were obtained from the HOCOMOCO, FactorBook, HOMER, and ENCODE themes provided by the “MotifDb” package. The effect of somatic mutations on TF binding motifs was measured using motifbreakR, which employs an efficient algorithm based on information content to distinguish between truly destructive and neutral variants. For each variant, motifbreakR first queries the PWM database to determine which motif the variant may belong to, and then returns metrics assessing the degree to which each allele of the variant disrupts (or enhances) the binding of the putative TF binding motif. Finally run the function motifbreakR(), using the parameters filter = TRUE, threshold = 1e-4, method = 'ic', bkg = c(A = 0.25, C = 0.25, G = 0.25, T = 0.25), and BPPARAM = BiocParallel:: SerialParam(). 59 , 60 Univariate cox regression analysis and multivariate cox regression analysis were implemented to identify mutation-driven enhancers that were independently associated with prognosis of breast cancer patients. According to clinical information on breast cancer cases in TCGA, the samples with survival time of 0 and missing were deleted. Univariate cox proportional hazards model was used to screen out enhancers that were significantly related to survival with p < 0.05 as the threshold. Then multivariate cox regression analysis was performed to determine the independent effect of risk mutation-driven enhancers for overall survival. The risk score of each sample was calculated according to the formula: risk score = ∑ i = 1 n c o e f i × e n h a n c e r i where coef i is the regression coefficient of enhancer i . The enhancer i represents the chromatin accessibility of enhancer i . The n represents the number of enhancers. Patients were divided into a high-risk group and a low-risk group according to the median of the risk score. Kaplan–Meier analysis was used to compare overall survival between a high-risk group and a low-risk group. To identify candidate small molecules that affect mutation-driven enhancers, we obtained ranked gene lists from drug-induced gene expression profiles based on the Broad Institute’s Connectivity Map. We calculated whether drug perturbations would significantly affect the target genes of enhancers by performing Gene Set Enrichment Analysis (GSEA) on each drug in the Connectivity Map. We considered a predicted drug of mutation-driven enhancers if the FDR-adjusted P-value was less than 5% (FDR<0.05). All statistical analyses in this research were performed using R software (version 4.2.0). Kaplan-Meier survival curves were plotted using “survival” in the R and the Log-rank test was used to compare OS between the two groups. Pearson test was used to compare the mRNA expression levels of mutant and non-mutant samples. The “Coxph” function in the “survminer” package was used for COX regression analysis. The “survminer” and “pROC” packages were used to calculate AUC for the ROC. All statistical tests were considered statistically significant at p < 0.05.