Comprehensive Analysis of Regulatory Factors and Immune-Associated Patterns to Decipher Common and BRCA1/2 Mutation-Type-Specific Critical Regulation in Breast Cancer

Introduction

Breast cancer (BC) is highly heterogeneous with potentially aggressive and complex biological features. Gene expression profiling has been used to analyze this complexity and provide new insights (Holm et al., 2012; Jonsson et al., 2012). Certain genetic mutations typically lead to distinct subtypes of tumors. Sporadic basal-like cancer has phenotypically similarity and relationships with BRCA (BC susceptibility gene)-associated BC (Richardson et al., 2006; Alimonti et al., 2010). BRCA1 and BRCA2 are tumor suppressor genes with critical roles in DNA damage repair; their germline or somatic mutations are closely related to carcinogenesis, including BC and high-grade serous cancer of the gynecological tract (particularly the fallopian tube, ovary, and peritoneum) (Foulkes and Shuen, 2013). BRCA1/2 mutations are the most common germline mutations and confer substantial lifetime risk of tumors, accounting for up to 40% of familial BC cases (Anastasiadi et al., 2017). Vidula et al. (2020) evaluated the characteristics of somatically acquired BRCA1/2 mutations in BC patients and found that 13.5% had somatic BRCA1/2 mutations (of which 4% were known germline pathogenic variants and the remainder novel variants) in metastatic BC. Germline variants of BRCA1/2 mutations in BC have been well described, whereas the influence of somatic variant profiles remains less clear. Owing to the introduction of PARP inhibitors and DNA damaging agents as therapies for germline BRCA1/2-mutant (BRCA1/2-MUT) advanced BC and germline/somatic BRCA-MUT ovarian cancer, the identification of non-familial BC cases with BRCA-like features is important, as these may represent a novel therapy opportunity (Tutt et al., 2018; Coleman et al., 2019; Vidula et al., 2020). Generally speaking, lacking the specific attention to the somatic mutations and effects in BRCA1/2 genes on the molecular carcinogenic mechanism and progression in BC.

The objectives of this study were to discover genomic expression and immune patterns closely associated with the carcinogenesis of BRCA mutations, based on BC sequencing from The Cancer Genome Atlas (TCGA) database, and to identify the similarities and specificities underlying disease pathogenesis of somatic BRCA-mutant and wildtype BC. This was expected to provide novel insights into the molecular mechanisms influenced by this disease and BRCA mutations, and to identify putative targets for preventive strategies to treat individuals with this tumor susceptibility gene.

Triple-negative BC (TNBC) is a particular type of BC that is characterized by its heterogeneity, early recurrence, and poor survival (Foulkes et al., 2010). TNBC has an important relationship with BRCA1/2-MUT BC in terms of pathology and immunophenotype (Sharma et al., 2014; Belli et al., 2019). Therefore, given the lack of effective targeted therapies, the poor prognosis of TNBC patients, and the close relationship between BRCA1/2-MUT BC and TNBC, we also aimed to further identify closely related genes between these two subtypes, in order to reveal underlying molecular relations and promising therapeutic targets.

Some studies on BC have obtained important insights using high-throughput data and weighted gene co-expression network analysis (WGCNA) (Deng et al., 2019). WGCNA is used to construct scale-free gene co-expression networks and to find modules of genes with similar expression patterns under specific circumstances, which provides insights into key genes and signaling networks that could play critical roles in the progression of diseases by linking clinical phenotypes. However, as traditional WGCNA extensively used in many diseases studies, requires datasets of similar biological significance and analysis of only one group (one subtype), it may not be a suitable technique for analyzing differences in biological functions among several subtypes of diseases. Consensus WGCNA (Consensus network analysis) is rarely used in practice after it was proposed, but it has revealed changes in pathway and molecular mechanisms dependencies due to differential biological subtypes and achieved biologically meaningful results, in some studies such as the male-female mouse liver tissues, and human-canine osteosarcoma (Langfelder and Horvath, 2007; Jin et al., 2019; McDonough et al., 2019). Therefore, to identify the differences between BRCA1/2-MUT BC and BRCA1/2 wild-type (BRCA1/2-WT) BC, we used consensus WGCNA, an improved analysis method which was conducted in our study using a slightly different R package code on the basis of the general principles of WGCNA, thereby perform the comparison of the effect of modules on associated clinical traits between mutant and WT cancers. Therefore, by Consensus network analysis, we have identified a series of important modules and genes within the modules, which have significantly universal/specific effects on the tumor stage, prognosis or metastasis indicators in BRCA1/2-MUT tumors, but some lacking obvious effects on wild-type tumors.

Immunophenotypes [estrogen receptor (ER); progesterone receptor (PR); human epidermal growth factor receptor 2 (HER2)] in BC, were of great significance for guiding clinical treatment and judging the treatment response, disease outcome, patient recurrence and prognosis. For instance, many studies demonstrated that Her-2 positive and TNBC subtypes had a poor prognosis and are more likely to relapse and metastasize early and frequently. The anatomic assessment of tumor size, regional lymph node involvement, and distant metastases (known as TNM), and metastasis indicators including number of positive lymph nodes, and cytokeratin, margin status as well, could provide the important predictors for the distinct clinical behavior and prognosis of BC (Hortobagyi et al., 2018). In our study, we focused on the above clinical traits of two BC subtypes, and performed the module-traits relationship analysis. This was followed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment, and ENCODE chromatin immunoprecipitation sequencing (ChIP-seq), to identify biological function and transcription factors (TFs) associated with modules. Highly connected genes within modules were identified and validated using TCGA and Gene Expression Omnibus (GEO) data, and survival analysis was performed, and shown to be BRCA1/2-associated critical regulatory genes. Furthermore, we identified BRCA1/2-MUT-specific immune patterns. Accordingly, our results would indicate common and specific patterns of regulation of tumor metastasis and the tumor immune microenvironment (TME) in BC, providing novel insights into the pathological process of this disease and the corresponding BRCA mutations.

Materials and Methods

Study Population

All data were downloaded from TCGA¹. BRCA1-MUT and BRCA2-MUT BC samples were selected using cBioPortal² (Gao et al., 2013) based on the presence of BRCA1/2 mutations and deep deletion copy number variation (CNV). We excluded BC samples with BRCA1/2 amplification to avoid confounding factors. BRCA1/2-WT samples were randomly selected from BC patients without BRCA mutations or CNV; we also selected para-carcinoma tissues from our last study as normal samples (Li et al., 2020). Comparison of BRCA1-MUT and BRCA2-MUT samples with normal samples identified 4,889 and 5,124 differentially expressed genes (DEGs), respectively, using the EdgeR package. Differentially expressed RNA-sequencing expression profiles were used as an input for WGCNA, selected from more than 19,000 original profiles by keeping only the DEGs.

In the subsequent validation of the results, we used TNBC data from TCGA for internal validation and GEO datasets GSE3744 (Richardson et al., 2006; Alimonti et al., 2010) and GSE25307 (Jonsson et al., 2012) for external validation. The research subjects used in our work and relevant clinical traits information are shown in Table 1.

TABLE 1

Table 1. Main subjects used in this study.

Weighted Gene Co-expression Network Analysis

Consensus WGCNA was conducted using the R package based on R version 3.5.1, to determine co-expressed genes between the two subtype groups from among the DEGs found in BRCA1/2-MUT BC. This was intended to identify key regulatory factors in the progression of BRCA1/2-MUT BC. We used fragments per kilobase of transcript per million fragments mapped (FPKM), or transcripts per million (TPM) normalization methods to standardize sequencing depth and gene length, in order to measure gene expression levels (Zhao et al., 2020). The MUT and WT expression profile data processed by FPKM were used for network construction. Hierarchical clustering of samples was used to identify outliers in each group (Figure 2). Based on the criterion of approximate scale-free topology, we chose the appropriate soft thresholding power for the function pickSoftThreshold; this is the power to which the co-expression similarity is raised to calculate adjacency (Supplementary Figure 3). We chose soft thresholding power values of 14 and 5, minimum module sizes 30 and 45, module detection sensitivity deepSplit 2, and cut height for merging of modules of 0.25, respectively, for BRCA1-MUT and BRCA2-MUT BC. Co-expression genes were assigned to modules via dynamic minimum tree-cutting arithmetic. Consensus network analysis was performed between MUT and WT patients (Supplementary Figure 4; refer to the research of Horvath et al., could get more detailed method introduction) (Langfelder and Horvath, 2007).

Identifying consensus modules significantly correlated with clinical traits in the MUT and WT subtypes, and with the same direction of correlation (positive or negative), could help to determine important modules correlated with clinical traits, including ER/PR and tumor stage. Upregulation and downregulation of each consensus module was determined based on the log2 fold change (log2FC) value: the mean gene expression in the mutant group was divided by the corresponding gene expression in normal tissue, and the log (base 2) of this value was taken. Module membership was determined by Pearson correlation between gene expression profiles and modules and used to rank module connectivity, in order to identify genes with high module membership in modules of interest, as well as those with high significance for traits.

Enrichment Analysis for Biological Functions and Pathways

Module biological function was determined using g:Profiler³ to analyze GO and signaling pathway enrichment. GO analysis was performed with respect to biological process, cellular component, and molecular function (Raudvere et al., 2019). Comprehensive biological pathway analysis was conducted using the online KEGG pathway, Reactomen, and WikiPathways resources. Enriched GO terms and pathways were considered to be significant according to the criterion of adjusted P < 0.05.

Identification of Associated Transcription Factors

Transcription factors for clinically significant modules (turquoise, purple, red, blue, turquoise, black, and brown) were identified using the iRegulon V.1.3 plugin in Cytoscape, based on 1120 ChIP-seq data from the ENCODE database, using the criterion of normalized enrichment score (NES) > 4.0 (Janky et al., 2014; McDonough et al., 2019).

Identification and Validation of Candidate Genes in the Cancer Genome Atlas and Gene Expression Omnibus Database

We identified highly connected genes within modules as candidate genes to analyze the changes in transcriptional expression levels in the BRCA1-MUT, BRCA2-MUT, and WT groups compared with normal samples. In the validation cohort, we performed DEG analysis between TNBC expression profiles and para-cancerous tissues to achieve internal validation of candidate genes; GSE3744 (Richardson et al., 2006; Alimonti et al., 2010) and GSE25307 (Jonsson et al., 2012) were downloaded for external validation, using the limma package and P < 0.05 was taken as the criterion for a significant difference. Furthermore, we analyzed and evaluated the mRNA expression levels of candidate genes using an R package combined with standardized TPM expression data for several subtype groups (sub-division into three groups with TNBC included or excluded) to reflect the potential effects of TNBC subtypes on the final results.

Survival Analysis of Candidate Genes for Breast Cancer

Candidate genes were evaluated with respect to their prognostic value in BC. We used Kaplan–Meier Plotter⁴, an online resource for survival analysis. As there were insufficient BRCA1/2-MUT BC cases for accurate survival analysis, the overall survival (OS) of all BC patients (N = 1879, from TCGA and GEO datasets) was analyzed using the Kaplan–Meier method, based on the classification of patients into high and low groups according to the mRNA expression of corresponding candidate genes. In the survival analysis, significance was defined as log-rank P < 0.05.

Correlation Analysis Between BRCA1/2-MUT BC and Infiltrating Immune Cells

The TIMER web server⁵ is a comprehensive resource for systematical analysis of immune infiltrates across diverse cancer types (Li et al., 2016, 2017). TIMER was used to perform correlation analysis between the expression levels of hub genes in our study and tumor purity, as well as levels of tumor-infiltrating immune cells including B cells, CD4 + T cells, CD8 + T cells, macrophages, neutrophils, and dendritic cells.

Correlations of Infiltrating Immune Cell Proportions, Tumor Mutational Burden, and Neoantigen Load With Hub Gene Expression

Based on bulk RNA-seq matrix information, bioinformatics methods using deconvolution-based analysis were used to score immune cells and the immune microenvironment, and to estimate immune cell infiltration. CIBERSORT enables characterization of cell composition and score estimation for 22 immune cell types in various tissues, based on gene expression profiles (Newman et al., 2015). We used the open-source CIBERSORT resource to compare infiltrating immune cell proportions in BC types with and without BRCA1/2 mutations (BRCA1-MUT, BRCA2-MUT, and BRCA1/2-WT), as well as in TNBC and normal tissues. Spearman correlation analysis was used to test the correlations between hub gene expression and immune infiltration in BRCA1/2-MUT BC. Relationships of gene expression with tumor mutational burden (TMB) and neoantigen load in all BC types were analyzed using the maftools R package with TCGA mutation data and the SangerBox database⁶, an online platform for TCGA patient data based on R (Hu et al., 2021; Zhang et al., 2021).

Statistical Analysis

Statistical analyses were performed with GraphPad Prism 7 and SPSS version 18.0. Differences among three or more groups were determined by analysis of variance or Kruskal–Wallis test, and further differences between groups were analyzed by the Dunnett method or student–Newman–Keuls test. Survival analysis was performed using the Kaplan–Meier method with log-rank test. Correlations between gene markers of tumor-infiltrating immune cells and hub genes were obtained as Spearman’s rho values. P < 0.05 was considered to indicate statistical significance.

Results

Weighted Gene Co-expression Network Analysis and Key Module Identification

Consensus network analysis identified 12 consensus modules (genes not assigned to any of the modules are colored gray), including eight upregulated modules and three downregulated modules, in the BRCA1-MUT cohort compared with the WT cohort (Figures 1A,B); and 12 consensus modules including nine upregulated modules and two downregulated modules in BRCA2-MUT BC samples (Figures 1C,D; a complete module-trait heatmap is provided in Supplementary Figure 4).

FIGURE 1

Figure 1. Modules associated with clinical traits in BRCA1/2-MUT and BRCA1/2-WT BC, and expression levels of corresponding module genes. (A,C) Identification of consensus modules for BRCA1/2-MUT and BRCA1/2-WT BC; complete figures are shown in Supplementary Figure 4. (B,D) Plot of gene expression levels for each consensus module of BRCA1-MUT and BRCA2-MUT BC to determine overall upregulation or downregulation of the module relative to normal tissue. N or M, the anatomic assessment of regional lymph node involvement or distant metastases. ER, status of estrogen receptor; PR, status of progesterone receptor; HER2, status of human epidermal growth factor receptor-2; lymph nodes, number of positive lymph nodes; ME, module eigengene. Negative ER, PR, or HER2 status was defined as “1” and positive status as “0.” Red frames, with significantly positive correlation in MUT; blue frames, lacking significant correlations in WT.

FIGURE 2

Figure 2. Network of enriched TFs (blue) and regulatory non-coding RNAs (yellow) with clinically significant modules (turquoise and purple modules in BRCA1-MUT BC; red, blue, turquoise, black, and brown modules in BRCA2-MUT). #ME represents modules from BRCA1-MUT samples. Red represents upregulated modules; green represents downregulated modules.

In the BRCA1-MUT group, ER/PR status was positively correlated with the turquoise module. Also, the turquoise and purple modules had significant correlations with tumor M stage, and the purple module had a significantly positive correlation with death, whereas no such correlations were found for WT samples (Figure 1A). But no obviously positive correlations were found for metastasis indicators in BRCA1-MUT samples (Supplementary Figure 4). These results suggested that genes in the turquoise and purple modules could explain the tumor metastasis mechanisms of BRCA1 mutations and thus influence clinical prognosis. In BRCA2-MUT samples, ER/PR status was correlated with eight modules (magenta, brown, greenyellow, pink, purple, blue, and black). N stage was positively correlated with the turquoise module. Metastasis indicators, including cytokeratin, positive lymph nodes, and margin, were significantly correlated with the turquoise, brown, red, greenyellow, purple, blue, and black modules, but there were no such significant correlations in the WT group (Figure 1C).

Enrichment Analysis for Biological Functions and Pathways

The turquoise and purple modules in the BRCA1-MUT group and the red, purple, blue, turquoise, black, greenyellow, and brown modules in the BRCA2-MUT group due to their correlations with clinical traits, were identified as clinically significant modules associated with carcinogenesis, progression, and metastasis of BRCA1/2-MUT BC. GO and pathway enrichment analyses were carried out for these modules (Table 2; complete results are shown in Supplementary Table 1).

TABLE 2

Table 2. Biological functions (from GO and pathway enrichment) in clinically significant modules.

In the BRCA1-MUT group, the turquoise module was enriched in cell cycle regulation, DNA replication, DNA damage response, and immune response; and the purple module was mainly enriched in processes related to the immune response, including defense response to virus and type I interferon signaling pathway, according to the GO enrichment analysis. Similarly, in the BRCA2-MUT group, modules were associated with the following categories: cell cycle, cell division, DNA replication, and DNA repair (blue module); cellular metabolic and catabolic processes (red module); cell adhesion and migration (red and brown modules); immune response (brown module); organ development (black module); and mitochondrial function (turquoise module). The purple and greenyellow modules showed no significant functional enrichment.

Identification of Module-Associated Regulators

In the comparison of BRCA1-MUT and BRCA1-WT BC samples, 17 TFs and five microRNAs (miRNAs) were found to be enriched in the turquoise and purple modules (Figure 2 and Supplementary Table 2). The strongest enrichment scores were obtained for TFs associated with the upregulated purple immune response (STAT2, STAT1, IRF4, and IRF1) and the upregulated turquoise cell cycle module (E2F7, E2F4, E2F1, FOXM1, and TFDP1), with all NES > 9. Based on their expression levels in BRCA1-MUT BC relative to normal samples, many TFs were also significantly upregulated in the BRCA1-MUT group, including in the purple module [STAT1 (log2FC = 1.62) and IRF1 (log2FC = 1.10)] and turquoise module [E2F7 (log2FC = 3.16), E2F1 (log2FC = 2.83), FOXM1 (log2FC = 4.13), E2F2 (log2FC = 3.03), and MYBL2 (log2FC = 4.20)], suggesting a positive feedback loop in the regulation of BRCA1-associated carcinogenesis.

In the BRCA2-MUT samples, 77 TFs and five miRNAs were enriched in the six clinically significant modules (Figure 2 and Supplementary Table 2). The important TFs with high NES values were as follows: E2F4 in the (blue) cell cycle/DNA repair module; SPIB and IRF9 in the (brown) immune response module; and STAT1 and EZR, among others, in the modules associated with cell migration, adhesion, or angiogenesis (red and brown). Also, many TFs were also significantly upregulated in the BRCA2-MUT group, including E2F1 (blue module, log2FC = 3.13), FOXM1 (blue module, log2FC = 3.82), STAT1 (red module, log2FC = 1.82), ZIC3 (turquoise module, log2FC = 2.59), and IRF9 (brown module, log2FC = 1.30), again suggesting the existence of a positive feedback loop.

In summary, the turquoise module in the BRCA1-MUT group and the blue module in the BRCA2-MUT group both showed specific associations with regulation of the cell cycle and DNA repair, and had correlations with miR-193b, miR-192, and miR-215 and enrichment of the E2F family. The purple module in the BRCA1-MUT group and the brown module in the BRCA2-MUT group were associated with regulation of important immune response functions and inflammatory pathways, in particular, members of the IRF family including IRF1, IRF4, and IRF9, which are associated with immune regulation.

Analysis of Clinically Significant Modules Including Functional Enrichment and Identification of Candidate Genes

The top 10 highly connected genes in each clinically significant module were screened for further analysis. As shown in Figures 3, 4, we determined the expression levels of these genes in the MUT, WT, and TNBC groups, as well in normal tissues for comparison. In several of the modules associated with immune response and cell cycle regulation, gene expression levels showed sequential changes from TNBC to mutant BC to WT BC. We found that some modules specifically associated with the clinical characteristics of BRCA2-MUT BC were related to tumor metabolism, cell migration, and other biological functions, and most candidate genes were also expressed in TNBC (Figure 4).

FIGURE 3

Figure 3. Biological function enrichment results for genes in clinically significant modules in BRCA1-MUT BC and gene network analysis. Biological function enrichment and corresponding gene information are shown for BRCA1-associated clinically significant modules. (A) Purple module, (B) Turquoise module. Left column: Top GO terms and KEGG pathways for genes from clinically significant modules. Middle column: Protein–protein interaction networks of top 10 most connected genes with highest degree in modules, constructed using STRING and Cytoscape. Right column: Heatmaps showing log2FC values for top 10 most connected genes. The ordinate is the gene name, and the abscissa (BRCA1-MUT, WT, and TNBC) represents different BC datasets from TCGA. Red represents log2FC > 0, blue represents log2FC < 0. Values are log2FC values of each gene in three sets, where FC is calculated as the mean gene expression within a group divided by the corresponding gene expression of normal tissues.

FIGURE 4

Figure 4. Biological function enrichment results for genes in clinically significant modules for BRCA2-MUT BC and gene network analysis. Biological function enrichment and corresponding gene information are shown for BRCA2-associated clinically significant modules. (A) Brown module, (B) blue module, (C) red module, (D) turquoise module, (E) black module. Left column: Top GO terms and KEGG pathways for genes from clinically significant modules. Middle column: Protein–protein interaction networks of top 10 most connected genes with highest degrees in modules, constructed using STRING and Cytoscape. Right column: Heatmaps showing log2FC values for top 10 most connected genes, similar to those shown in Figure 3.

Immune Response

In both the BRCA1-MUT and BRCA2-MUT cohorts, we found upregulated modules significantly related to immune response (the purple and brown modules, respectively), indicating that genes in these modules may have specific effects on progression and metastasis of BRCA1/2-MUT BC but no significant impact on BRCA1/2-WT BC. In the BRCA1-MUT group, the purple module was enriched in processes involved in immune response, including defense response to virus and regulation of cytokine production. On the other side, the pathway enrichment analysis results suggested that genes in this module were mainly related to the RIG-I-like receptor signaling pathway, JAK–STAT signaling pathway, and NOD-like receptor signaling pathway (Figure 3A). The top 10 genes ranked by degree were identified as candidate genes. These candidate genes included ISG15 and STAT1 (Figure 3A). The top seven genes had the highest node degrees among the candidate genes (i.e., node degree 31). In the BRCA2-MUT group, the brown module was enriched in similar GO biological process and pathways (Figure 4A), and CTLA4 had the highest node degree of 100.

Cell Cycle/DNA Repair

In the BRCA1-MUT group, the turquoise module was enriched in processes including cell cycle, cell cycle phase transition, and DNA repair (Figure 3B). Among the top 10 genes, which included BUB1 and CCNB1, BUB1 had the highest node degree (i.e., 150). In the BRCA2-MUT group, the blue module was enriched in similar GO biological processes and pathways (Figure 4B). Among the top 10 hub genes, which included CDK1, CCNB1, CCNA2, and BUB1, CDK1 had the highest node degree of 141.

Cell Migration/Adhesion, Metabolic Processes, and Organ Development

Notably, the red and brown modules (Figures 4C,D) were found to specifically influence the progression and metastasis of BRCA2-MUT BC, via their important roles in the regulation of cell migration and adhesion. No such correlations were observed in WT BC. The red module was also enriched in metabolic processes including regulation of lipolysis, response to insulin, gluconeogenesis, inflammatory response, and regulation of cell migration. On the other side, pathway enrichment analysis suggested that the genes in the red module were mainly related to the PPAR signaling, AMPK signaling, and estrogen signaling pathways (Figure 4C). The turquoise module was mainly associated with mitochondrial translation functions (Figure 4D), and the genes in the black module appeared to have roles in anatomical structure morphogenesis and organ development (Figure 4E).

Validation and Identification of Hub Genes

By internal validation using TCGA data, we demonstrated that most candidate genes were also significantly expressed in TNBC (Figures 3, 4). This indicated the accuracy of our results; the important roles of these genes in the progression of BC, especially in BRCA1/2-MUT BC; and the relationships between TNBC and BRCA1/2 mutations at the molecular level.

Moreover, we used GEO datasets GSE3744 and GSE25307 for further validation of these hub genes. As shown in Figure 5, there were many genes that were not found in the GEO datasets, and most other genes had differential expression relative to normal tissues. In BRCA1-MUT BC, genes including BUB1, CCNB1, BUB1B, CCNA2, and TTK showed consistently upregulated expression in at least one of the validation datasets. In BRCA2-MUT BC, ITGAX, CCNB1, CCNA2, BUB1, and BUB1B also showed upregulated expression, whereas LIPE and FABP4 were downregulated.

FIGURE 5

Figure 5. Heatmap of log2FC values for external validation of hub genes using GEO datasets and transcriptional expression levels of validated hub genes. (A,B) Heatmaps of log2FC values for hub genes in BRCA1-MUT (A) and BRCA2-MUT (B) BC cohorts for GEO external validation. The ordinate is the gene name and the abscissa represents different groups in GEO datasets. Values are log2FC values of each gene in three sets, where FC is calculated as the mean gene expression in the relevant group divided by the corresponding gene expression in normal tissue. “NA” represents GEO datasets that did not contain these genes. (C–E) Expression levels of validated hub genes in BRCA1-MUT (C,D) and BRCA2-MUT (C,E) BC. To show the potential effects of TNBC subtypes on BC, the hub gene expression was analyzed in multiple groups: normal (N = 112), BRCA1-MUT BC (N = 36), BRCA1-MUT BC excluding TNBC (N = 26), BRCA2-MUT BC (N = 49), BRCA2-MUT BC excluding TNBC (N = 42), TNBC (N = 116), WT (N = 117), and WT excluding TNBC (N = 103). **P < 0.01, ****P < 0.0001.

Given the different immune profile of the TNBC subtype compared with other BC subtypes, and the associations between BRCA1/2 mutations and TNBC, we selected non-TNBC tissues with BRCA1/2 mutations and non-TNBC BRCA1/2-WT tissues (that is, MUT and WT samples excluding all TNBC subtypes) for comparison. As shown in Figures 5C,D, BRCA1-associated hub genes BUB1, BUB1B, CCNB1, CDC20, IFI44, ISG15, KIF11, STAT1, RSAD2, TTK, and NCAPG showed upregulated gene expression. Similarly, BRCA2-associated hub genes BUB1, BUB1B, CCNB1, FABP4, and LIPE showed differential mRNA expression levels in BRCA2-MUT and non-TNBC BRCA2-MUT tissues, suggesting that TNBC-subtype effects did not completely prevent the identification of many hub genes. It is worth noting that in the TNBC subtype, some oncogenes, including BUB1, were further upregulated at the transcriptional level, whereas tumor suppressor genes FABP4 and LIPE were downregulated, compared with normal tissues, BRCA1/2-WT, and BRCA1/2 mutants.

In this study, we found that BUB1, BUB1B, CCNB1, CDC20, ISG15, KIF11, NCAPG, and TTK were not only significantly expressed in corresponding tumor tissues with BRCA somatic mutations compared with WT BC and normal tissue but also showed significant prognostic value for BC (Figures 6A–I). Therefore, these validated genes are potential key hub genes and critical regulatory factors in BRCA1/2-MUT BC, which could represent new prognostic or diagnostic biomarkers and provide novel insights into the mechanisms of BRCA1/2 mutations in cancer development.

FIGURE 6

Figure 6. Prognostic significance of validated hub genes in BC. (A–H) Validated hub genes with significant prognostic value for OS in BC are shown: (A) BUB1, (B) BUB1B, (C) CCNB1, (D) CDC20, (E) ISG15, (F) KIF11, (G) NCAPG, and (H) TTK. (I) Prognostic value of these eight hub genes for 5-year OS in BC.

Correlation Analysis Between Hub Genes of BRCA1/2-MUT BC and Infiltrating Immune Cells, Tumor Mutational Burden, and Neoantigen Load

We found that the validated hub genes could influence metastasis of BRCA1/2-MUT BC through mechanisms including regulation of immune response, cell cycle, cell migration, and metabolic processes. Many recent studies have demonstrated that tumor-infiltrating immunocytes (TICs) can regulate and participate in the TME, regulate tumor progression, and influence lymph node metastasis and prognosis (Azimi et al., 2012; Yusuf et al., 2020). Thus, we used the TIMER web server and CIBERSORT to assess immune infiltration and its correlations with BRCA1/2-associated hub genes.

We noted that hub gene expression was positively correlated with infiltrating immune cells, including T cells, neutrophils, and dendritic cells, in almost all BC subtypes (basal-like, luminal, and HER2-positive; Figure 7 and Supplementary Figure 5), of which the basal-like and luminal subtypes showed significant associations with tumor purity. ISG15 expression was negatively correlated with tumor purity in basal-like BC, as well as with levels of infiltrating T cells, neutrophils, and dendritic cells; its correlation with tumor purity in the luminal subtype was undefined (Figure 7A). However, BUB1, CCNB1, BUB1B, KIF11, CDC20, TTK, and NCAPG (Figures 7B–H), which are closely related to regulation of the cell cycle and DNA repair, showed strong positive correlations with tumor purity in luminal-type BC.

FIGURE 7

Figure 7. Correlation of each hub gene’s expression with immune infiltration levels in BC. Complete correlation analysis results for BC subtypes are given in Supplementary Figure 5. Correlation of ISG15 (A), BUB1 (B), CCNB1 (C), BUB1B (D), KIF11 (E), CDC20 (F), TTK (G), and NCAPG (H) expression with TICs.

Furthermore, to explore the effects of hub genes on TICs, we analyzed the relationships between their expression levels and various specific gene markers of immune cells, including innate immune cells and adaptive immune cells, after adjustment for tumor purity (Figure 8A and Supplementary Table 3). As shown in Figure 8A, we first determined the correlation coefficients between the expression levels of hub genes (ISG15, BUB1, CCNB1, BUB1B, KIF11, CDC20, TTK, and NCAPG) and TIC markers. Positive correlations with hub gene expression were observed for markers of Th1 cells including T-bet, STAT4, and IFN-r; markers of T regulatory cells (Tregs) including FOXP3 and CCR8; and markers of Tfh cells including IL21, CD278, and CXCL13. Markers of natural killer cells (NK cells), tumor-associated macrophages, and M1/M2 macrophages also showed possible positive correlations with expression of hub genes. We observed significant positive correlations with markers of T cell exhaustion, including FOXP3, PD-1, CTLA4, LAG3, TIM-3, and GZMB. We also found that BRCA1/2-specific infiltration patterns displayed a consistent lack of correlation with Th2 and Th17 cell markers.

FIGURE 8

Figure 8. Analysis and identification of BRCA1/2-MUT-specific immune-associated patterns and their relationships with TMB and neoantigen load in BC. (A) Correlation analysis between related hub genes (ISG15, BUB1, CCNB1, BUB1B, KIF11, CDC20, TTK, and NCAPG) and immune cell markers in TIMER. Expression scatter plots for pairs of genes adjusted by tumor purity in BC (1093 cases), together with Spearman’s rho value and estimated statistical significance, reflecting the correlations between hub genes and immune checkpoints/gene markers of immune cells. (B) Correlations between hub gene expression and infiltrating immune cell proportions in BRCA1-MUT and BRCA2-MUT BC. (C) Infiltrating immune cell proportions (with obvious correlations with hub gene expression in (B) in BRCA1/2-MUT BC, TNBC, BRCA1/2-WT BC, and para-cancerous normal tissues, using CIBERSORT package. (D,E) Comparison of immune cell proportion between stages I, and stage II + III and stage IV in BRCA1/2-MUT BC patients, and para-cancerous normal tissues. (F) Correlations between hub gene expression and TMB in BC. (G) Correlations between hub gene expression and neoantigen load in BC. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Comparison of cell proportions between BRCA1/2-MUT BC and normal tissues, and their correlations with hub genes (Figures 8B,C), also partially supported the idea that BRCA1/2 mutations could increase proportions of infiltrating immune CD4+ T cells, for instance, Th1, Tfh, and Treg cells, in BC. Specifically, the expression levels of some hub genes were positively correlated with levels of infiltrating CD4+ T cells, Tfh cells, Tregs, and M1 macrophages (Figure 8B). We also observed significantly higher levels of immune infiltration in BRCA1/2-MUT BC, especially for CD4+ T cells, Tfh cells, Tregs, and M1 macrophages, compared with normal tissue or BRCA1/2-WT BC (Figure 8C). This was followed by the further analysis of the above infiltrating immune cells proportions in different clinical stage of tumors. The results show that Tfh cells and Tregs have a higher proportion in patients with high tumor stage than in stage I, but decreased level trend of M0, M1 macrophages infiltration, compared with that of stage I tumors (Figures 8D,E), suggesting that with the progression of BRCA1/2-mutant BC, the proportions of infiltrating immune cells showed an increased trend, especially Tfh and Tregs.

The high TMB and neoantigen load in tumors favor the infiltration of immune effector cells. Antitumor immunotherapy responses are strong in these patients relative to others (Picard et al., 2020). In many cancer types, higher TMB may be associated with favorable response to anti-PD-1/PD-L1 immunotherapy (Wang and Li, 2019). In our study, higher expression levels of all hub genes were positively associated with TMB in BC patients (Figure 8F). ISG15, CDC20, TTK, and NCAPG expression showed positive correlations with neoantigen load in TCGA BC patients (Figure 8G). These results suggest that the hub genes, especially ISG15, CDC20, TTK, and NCAPG, could predict the efficacy of immunotherapy in BC.

Therefore, these complementary results demonstrate strong positive correlations of BRCA1/2-associated hub genes with T cell infiltration and T cell exhaustion, especially that involving T cell subtype markers such as Th1 and Tfh cells and Tregs. These findings suggest that BRCA1/2 mutation of breast tumors could recruit immune cells in the TME, resulting in the recruitment of Th1 cells, Tfh cells, and Tregs, and T cell exhaustion.

Discussion

Breast cancer remains the most common cancer affecting female patients and a major cause of cancer-related deaths among women worldwide (DeSantis et al., 2014). There are ongoing efforts to assess the mutational profiles of BC. Somatic mutations in TP53 and PIK3CA occur at rates of about 25% (Azim and Partridge, 2014; Weisman et al., 2016); on the other hand, there are numerous novel subtype-associated genetic mutations, including specific mutations of GATA3, PIK3CA, and MAP3K1 in the luminal A subtype (Cancer Genome Atlas Network, 2012). However, genetic evidence regarding mutation profiles and their implications for female patients remain to be clarified. There were few related studies about explaining the influence of specific tumor somatic mutation types on the molecular portraits in BC. Owing to the need to explore prognostic or predictive markers and therapeutic targets for diseases, the prevailing research strategies include use of large-scale genome sequencing data combined with bioinformatic analyses (Wu et al., 2019). The transcriptome landscapes of specific subtypes of tumors represent a promising area of study; given the trend of tumor precision medicine, this has become an important direction of bioinformatics research.

In this study, we constructed a consensus co-expression gene network of DEGs influenced by BRCA1/2 mutations and associated TFs, based on results from our previous study (Li et al., 2020). The aim of the present study was to explore the regulatory factors and immune-associated patterns that have important common and specific influences and are involved in critical regulation of BRCA1/2 mutations. This was expected to improve our understanding of the molecular mechanisms of BC and provide insights for use in individualized treatment. Consensus WGCNA has rarely been used in research on this related subject. Consensus WGCNA enabled us to focus on the differences between the several subtypes with respect to transcriptome characteristics associated with clinical traits, in order to identify the similarities and differences between several subtypes of BC associated with tumor somatic mutations or prognosis/metastasis types.

The two disease subtypes (genetic mutation and wild type) were biologically very similar, but significant differences were observed, including multi-gene and multi-signal differences, as demonstrated by consensus network analysis in our study. Consensus network analysis identified the specific and common carcinogenic mechanisms involving BRCA1/2 mutations in BC with relation to differential clinical traits. For example, we found that two modules associated with immune response and cell cycle regulation were significantly correlated with clinical traits (tumor M stage and poor prognosis) in BRCA1-MUT BC, but without significant influence in WT BC, suggesting that the modules and associated genes influenced by BRCA1 mutations could regulate and cause dysfunction of the cell cycle and immune response, and thus influence the metastasis and prognosis of BRCA1-MUT BC. In BRCA2-MUT BC, the critical co-expression modules were significantly positively correlated with N stage, cytokeratin, number of positive lymph nodes, and margin indicators but lacked significant correlations in WT BC. Therefore, we identified these modules as BRCA2-MUT-type-specific co-expressed modules. As number of positive lymph nodes, cytokeratin, and margin indicators are important clinical prognostic factors, as demonstrated by many studies, these results suggested that BRCA1/2 could influence the genes in these modules and thus affect clinical outcomes (Engel et al., 2019). In addition, specific BRCA2-associated modules showed obvious effects on regulation of the cell cycle, DNA repair, cell adhesion/migration, immune response, and cellular metabolic processes via enriched pathways including p53 signaling, PPAR signaling, the JAK–STAT pathway, and the PD-1 checkpoint pathway.

Our results were consistent with some of the key functions of BRCA proteins, including cell cycle transition and the DNA damage response, as demonstrated by previous studies. For instance, BRCA1 mutations have been shown to trigger deregulation of genes involved in the G2/M cell cycle transition, whereas BRCA2 mutations could cause deregulation of genes involved in the G1/S cell cycle checkpoint (Bouwman and Jonkers, 2012; Vaclova et al., 2015; Tanaka et al., 2018). However, there was a lack of research evidence for other enriched biological processes and signaling pathways identified in our study, such as immune response and the JAK–STAT pathway.

Cytokeratin and tumor margin indicators are essential for the development and metastasis of various cancers and are regarded as important predictors of distant metastasis risk and poor prognosis (Lee et al., 2019; Li et al., 2019). A recent study found that BRCA1/2 carriers frequently experienced lung, distant lymph node, and central nervous system involvement (Song et al., 2020). In our study, the co-expressed gene network displayed stronger relationships with positive lymph nodes, cytokeratin, and margin indicators, compared with WT BC, suggesting that hub genes influenced by BRCA proteins could more effectively promote metastasis. Regarding the possible molecular mechanisms underlying this effect, our results suggest that BRCA1 mutations could influence cell cycle transition and tumor immunomodulation through the p53 or JAK–STAT pathway, whereas BRCA2 mutations could influence cell cycle transition, tumor immunomodulation cell migration, and energy metabolic processes via PPAR signaling, AMPK signaling, p53 signaling, the JAK–STAT pathway, the NF-κB pathway, and the regulation of cell adhesion molecule interactions. These effects could be related to metastasis risk and poor prognosis in BC patients.

Furthermore, regulators cross-linking these modules include the transcription factors STAT1, IRF1, E2F7, E2F1, E2F2, FOXM1, and MYBL2 which were significantly upregulated in BRCA1-MUT BC, and E2F1, FOXM1, STAT1, ZIC3, and IRF9 upregulated in BRCA2-MUT BC. These results suggested the existence of a positive feedback loop involving regulatory factors but also indicate potential therapeutic targets.

The ultimate goal of precision medicine is to identify specific molecular alterations that permit application of effective targeted drugs (Belli et al., 2019). Many studies have focused on the relationships between BRCA mutations and TNBC, and BRCA1/2 mutations are expected to provide potential therapeutic targets for TNBC. In this respect, although PARP inhibitors have been approved for the treatment of BRCA-MUT and invasive BC in HER2-negative patients, the emergence and development of drug resistance is inevitable. Therefore, more research on therapeutic targets is needed with respect to the precision treatment of TNBC and BRCA-MUT BC. The genes influenced by BRCA mutations will provide novel insights into the pathological processes underlying disease mechanisms in BRCA-MUT BC and guide therapeutic strategies. Notably, the highly connected genes identified in our study, denoted hub genes, showed sequential changes in expression levels from TNBC to BRCA1/2-MUT BC to BRCA1/2-WT BC, suggesting close relationships between BRCA1/2 mutations and TNBC at a molecular level. Furthermore, by validation using GEO datasets and survival analysis, we identified hub genes BUB1, CCNB1, BUB1B, ISG15, KIF11, CDC20, TTK, and NCAPG, which act as important regulatory genes in the progression of BRCA1/2-MUT BC.

In fact, some studies have indicated that BUB1 (Takagi et al., 2013), KIF11 (Pei et al., 2017), TTK (Sugimoto et al., 2017) and so on (Yuan et al., 2006) might be involved in tumorigenesis and become new candidate biomarkers for BC or TNBC treatment. Recent studies found that ISG15 could become markers predicting the metastasis in TNBC (Rojas et al., 2019). NCAPG was a newly found markers associated with the prognosis of BC (Xiao et al., 2020), but its mechanism and more evidence needed to be further reported. And the important functions of these genes in BRCA1/2-mutant BC were still lack of relevant studies. In our study, the upregulated expression and prognostic value of these hub genes in BRCA1/2-mutant BC suggested that its roles in this type of breast cancer were equally noteworthy. More importantly, we found these hub genes were closely associated with higher levels of immune infiltration in BRCA1/2-mutant BC, especially CD4+ T cells, Tfh cells, Tregs, and M1 macrophages and T cell exhaustion. ISG15, CDC20, TTK, and NCAPG expression showed positive correlations with tumor mutational burden and neoantigen load in BC patients. Our study indicated the importance of further exploration of these hub genes with potential to recognize immune infiltrated and predict the efficacy of immunotherapy.

Characterizing the mechanistic relationships of cancer with inflammation and the immune microenvironment can be conducive to identifying novel individualized biomarkers or therapeutic targets. The TME is essential in the development and progression of cancer. Immune response and inflammation are classic and prevalent examples of functions of the TME. Moreover, evasion of immune destruction and tumor-promoting inflammation are essential hallmarks of cancer (Hanahan and Weinberg, 2011; Iyengar et al., 2016; Zhao et al., 2019). The relationships among immune response, inflammation, and BC are complex and important. The roles of TILs within the immune microenvironment have received increasing attention because of their important regulatory effects on the progression, metastasis, and prognosis of tumors. Our functional enrichment results demonstrated that multiple genes affected by BRCA1/2 mutations showed significant enrichment in immune response and lymphocyte activation in BC. Some other studies reported that TNBC was the most immunogenic subtype and attracted TILs owing to its genomic instability and higher mutation rate (Oshi et al., 2020). However, the global characteristics of the immune response and inflammation in BC with BRCA1/2 mutations remained unknown.

Our study indicated a meaningful positive correlation of BRCA1/2-associated hub genes with infiltration of T cells (Th1, Tfh, and Treg cells) and T cell exhaustion. Engel et al. (2014) reported that infiltration of Tregs was significantly higher in TNBC or in BC with mutation of BRCA1. A recent study indicated that Treg frequencies increased with nodal invasion, and that Treg-mediated local immunosuppression might influence tumor-draining lymph node invasion by metastatic cells and lead to poor prognosis (Núñez et al., 2020). There is also evidence that Tregs may negatively regulate responses of CD8+ T cells and natural killer cells to tumor cells, as well as leading to severe neovascularization and metastasis, together with the effects of Tfh and B cells, as a pro-tumor immune response (Marshall et al., 2016). Many studies have shown that Th1 cells play a critical part in anti-tumor immunity through inducing and directing effector cellular immune responses as well as inflammatory responses, together with their cytokines (especially IFNγ), in several solid tumor types (Faghih et al., 2014). It is discovered that the original Tfh cells in BC promote the accumulation of Tregs in the tumor, inhibit the anti-tumor immunity, and ultimately promote the development of BC (Gu-Trantien et al., 2017). Notably, our study also revealed a series of further and detailed immune infiltration information influenced by BRCA1/2 mutation type in BC. Our findings implied that BRCA1/2 mutations could regulate multiple genes and pathways as described above, resulting in specific immune patterns, with effects on the TME including activation and infiltration of Tregs, Th1 cells, and Tfh cells, as well as T cell exhaustion. Therefore, our study demonstrates a potential mechanism by which BRCA1/2 could regulate T cell function and influence infiltrating immune cell recruitment and functional suppression in the TME. Furthermore, our findings indicated that during the development of BRCA1/2-mutant BC, Tfh cells and Tregs played a tumor-promoting role, and their proportions showed increased levels, suggesting the possibility that immune escape associated with the regulation of Tregs was aggravated, and the tumor’s ability to fight the immune response might further strengthened. However, we also observed the proportions of NK cells, M1 macrophages with the potential to exert an antitumor effect showed dynamic changes, which indicate the opposition between antitumor immune cells and tumor-promoting immune cells at a proportional level, and the escalating immune confrontation in the progression of BRCA1/2-mutant BC.

One important limitation of the study was that the BRCA1/2-MUT BC data did not distinguish between somatic mutations and germline mutations of BRCA genes. Therefore, we were unable to confirm whether the findings would be applicable to hereditary BC patients. Although duplicate verifications using TCGA and GEO data were performed to identify more credible key genes, further molecular biological experiments are required to support this investigation and confirm our results.

Conclusion

Our data demonstrated enrichment of several pathways associated with BRCA mutations; these results were partly consistent with current knowledge of the functions of BRCA. We identified a correlation gene expression network and associated TFs in BRCA1/2-MUT BC to provide comprehensive and novel insights into the pathological mechanisms involved in BRCA1/2-MUT BC and the close relations between BRCA1/2 mutations and TNBC at a molecular level. Our study identified several candidate regulatory factors and hub genes as promising targets for intervention, including members of the E2F TF family and IRF family. Moreover, BUB1, CCNB1, BUB1B, ISG15, KIF11, CDC20, TTK, and NCAPG showed prognostic value and potential roles in assessment of the efficacy of BC immunotherapy. These results indicate the close relationships between BRCA1/2 mutations and TNBC at the molecular level. Furthermore, we identified immune-associated patterns in BRCA1/2-MUT BC, which manifested as the recruitment of T cells, including Th1 and Tfh cells, Treg infiltration, and T cell exhaustion; these patterns would explain the critical influence on tumor progression and metastasis of BRCA1/2 mutations in BC via common and specific TME regulation patterns.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

YL and QH conceived and designed the experiments. WD collated and checked the data. YL and XZ performed the consensus WGCNA network. PZ and TZ performed the function enrichment analysis. LM, MQ, and PZ helped to perform the analysis about tumor-infiltrated immune cells. XM explained and collected the content about BRCA1/2 mutations. YL was responsible for bioinformatics and bio-statistics analysis, and wrote the manuscript. QH revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Natural Science Basic Research Program of Shaanxi Province (2020JZ-40).

Conflict of Interest

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

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.750897/full#supplementary-material

Footnotes

1. ^ https://portal.gdc.cancer.gov/
2. ^ http://www.cbioportal.org/index.do
3. ^ https://biit.cs.ut.ee/gprofiler
4. ^ http://kmplot.com/analysis/
5. ^ https://cistrome.shinyapps.io/timer/
6. ^ http://www.sangerbox.com/

References

Alimonti, A., Carracedo, A., Clohessy, J. G., Trotman, L. C., Nardella, C., Egia, A., et al. (2010). Subtle variations in Pten dose determine cancer susceptibility. Nat. Genet. 42, 454–458. doi: 10.1038/ng.556

PubMed Abstract | CrossRef Full Text | Google Scholar

Anastasiadi, Z., Lianos, G. D., Ignatiadou, E., Harissis, H. V., and Mitsis, M. (2017). Breast cancer in young women: an overview. Updates Surg. 69, 313–317. doi: 10.1007/s13304-017-0424-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Azim, H. A. Jr., and Partridge, A. H. (2014). Biology of breast cancer in young women. Breast Cancer Res. 16:427. doi: 10.1186/s13058-014-0427-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Azimi, F., Scolyer, R. A., Rumcheva, P., Moncrieff, M., Murali, R., McCarthy, S. W., et al. (2012). Tumor-infiltrating lymphocyte grade is an independent predictor of sentinel lymph node status and survival in patients with cutaneous melanoma. J. Clin. Oncol. 30, 2678–2683. doi: 10.1200/jco.2011.37.8539

PubMed Abstract | CrossRef Full Text | Google Scholar

Belli, C., Duso, B. A., Ferraro, E., and Curigliano, G. (2019). Homologous recombination deficiency in triple negative breast cancer. Breast 45, 15–21. doi: 10.1016/j.breast.2019.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouwman, P., and Jonkers, J. (2012). The effects of deregulated DNA damage signalling on cancer chemotherapy response and resistance. Nat. Rev. Cancer 12, 587–598. doi: 10.1038/nrc3342

PubMed Abstract | CrossRef Full Text | Google Scholar

Cancer Genome Atlas Network (2012). Comprehensive molecular portraits of human breast tumours. Nature 490, 61–70. doi: 10.1038/nature11412

PubMed Abstract | CrossRef Full Text | Google Scholar

Coleman, R. L., Fleming, G. F., Brady, M. F., Swisher, E. M., Steffensen, K. D., Friedlander, M., et al. (2019). Veliparib with first-line chemotherapy and as maintenance therapy in ovarian cancer. N. Engl. J. Med. 381, 2403–2415. doi: 10.1056/NEJMoa1909707

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, J. L., Xu, Y. H., and Wang, G. (2019). Identification of potential crucial genes and key pathways in breast cancer using bioinformatic analysis. Front. Genet. 10:695. doi: 10.3389/fgene.2019.00695

PubMed Abstract | CrossRef Full Text | Google Scholar

DeSantis, C., Ma, J., Bryan, L., and Jemal, A. (2014). Breast cancer statistics, 2013. CA Cancer J. Clin. 64, 52–62. doi: 10.3322/caac.21203

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, J. B., Honig, A., Kapp, M., Hahne, J. C., Meyer, S. R., Dietl, J., et al. (2014). Mechanisms of tumor immune escape in triple-negative breast cancers (TNBC) with and without mutated BRCA 1. Arch. Gynecol. Obstet. 289, 141–147. doi: 10.1007/s00404-013-2922-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, J., Weichert, W., Jung, A., Emeny, R., and Holzel, D. (2019). Lymph node infiltration, parallel metastasis and treatment success in breast cancer. Breast 48, 1–6. doi: 10.1016/j.breast.2019.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Faghih, Z., Erfani, N., Haghshenas, M. R., Safaei, A., Talei, A. R., and Ghaderi, A. (2014). Immune profiles of CD4+ lymphocyte subsets in breast cancer tumor draining lymph nodes. Immunol. Lett. 158, 57–65. doi: 10.1016/j.imlet.2013.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Foulkes, W. D., and Shuen, A. Y. (2013). In brief: BRCA1 and BRCA2. J. Pathol. 230, 347–349. doi: 10.1002/path.4205

PubMed Abstract | CrossRef Full Text | Google Scholar

Foulkes, W. D., Smith, I. E., and Reis-Filho, J. S. (2010). Triple-negative breast cancer. N. Engl. J. Med. 363, 1938–1948. doi: 10.1056/NEJMra1001389

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Aksoy, B. A., Dogrusoz, U., Dresdner, G., Gross, B., Sumer, S. O., et al. (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci. Signal. 6:l1. doi: 10.1126/scisignal.2004088

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu-Trantien, C., Migliori, E., Buisseret, L., de Wind, A., Brohée, S., Garaud, S., et al. (2017). CXCL13-producing TFH cells link immune suppression and adaptive memory in human breast cancer. JCI Insight 2:e91487. doi: 10.1172/jci.insight.91487

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of cancer: the next generation. Cell 144, 646–674. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Holm, K., Staaf, J., Jonsson, G., Vallon-Christersson, J., Gunnarsson, H., Arason, A., et al. (2012). Characterisation of amplification patterns and target genes at chromosome 11q13 in CCND1-amplified sporadic and familial breast tumours. Breast Cancer Res. Treat. 133, 583–594. doi: 10.1007/s10549-011-1817-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hortobagyi, G. N., Edge, S. B., and Giuliano, A. (2018). New and important changes in the TNM staging system for breast cancer. Am. Soc. Clin. Oncol. Educ. Book 38, 457–467. doi: 10.1200/edbk_201313

CrossRef Full Text | Google Scholar

Hu, J., Qiu, D., Yu, A., Hu, J., Deng, H., Li, H., et al. (2021). YTHDF1 is a potential pan-cancer biomarker for prognosis and immunotherapy. Front. Oncol. 11:607224. doi: 10.3389/fonc.2021.607224

PubMed Abstract | CrossRef Full Text | Google Scholar

Iyengar, N. M., Gucalp, A., Dannenberg, A. J., and Hudis, C. A. (2016). Obesity and cancer mechanisms: tumor microenvironment and inflammation. J. Clin. Oncol. 34, 4270–4276. doi: 10.1200/jco.2016.67.4283

PubMed Abstract | CrossRef Full Text | Google Scholar

Janky, R., Verfaillie, A., Imrichova, H., Van de Sande, B., Standaert, L., Christiaens, V., et al. (2014). iRegulon: from a gene list to a gene regulatory network using large motif and track collections. PLoS Comput. Biol. 10:e1003731. doi: 10.1371/journal.pcbi.1003731

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Z., Liu, S., Zhu, P., Tang, M., Wang, Y., Tian, Y., et al. (2019). Cross-species gene expression analysis reveals gene modules implicated in human osteosarcoma. Front. Genet. 10:697. doi: 10.3389/fgene.2019.00697

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonsson, G., Staaf, J., Vallon-Christersson, J., Ringner, M., Gruvberger-Saal, S. K., Saal, L. H., et al. (2012). The retinoblastoma gene undergoes rearrangements in BRCA1-deficient basal-like breast cancer. Cancer Res. 72, 4028–4036. doi: 10.1158/0008-5472.can-12-0097

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2007). Eigengene networks for studying the relationships between co-expression modules. BMC Syst. Biol. 1:54. doi: 10.1186/1752-0509-1-54

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, W. K., Lee, J., Kim, H., Lee, S. G., Choi, S. H., Jeong, S., et al. (2019). Peripheral location and infiltrative margin predict invasive features of papillary thyroid microcarcinoma. Eur. J. Endocrinol. 181, 139–149. doi: 10.1530/eje-18-1025

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., Severson, E., Pignon, J. C., Zhao, H., Li, T., Novak, J., et al. (2016). Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 17:174. doi: 10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Chen, Q., Deng, Z., Chen, X., Liu, H., Tao, Y., et al. (2019). KRT17 confers paclitaxel-induced resistance and migration to cervical cancer cells. Life Sci. 224, 255–262. doi: 10.1016/j.lfs.2019.03.065

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 77, e108–e110. doi: 10.1158/0008-5472.can-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Zhou, X., Liu, J., Yin, Y., Yuan, X., Yang, R., et al. (2020). Differentially expressed genes and key molecules of BRCA1/2-mutant breast cancer: evidence from bioinformatics analyses. PeerJ 8:e8403. doi: 10.7717/peerj.8403

PubMed Abstract | CrossRef Full Text | Google Scholar

Marshall, E. A., Ng, K. W., Kung, S. H., Conway, E. M., Martinez, V. D., Halvorsen, E. C., et al. (2016). Emerging roles of T helper 17 and regulatory T cells in lung cancer progression and metastasis. Mol. Cancer 15:67. doi: 10.1186/s12943-016-0551-1

PubMed Abstract | CrossRef Full Text | Google Scholar

McDonough, J. E., Kaminski, N., Thienpont, B., Hogg, J. C., Vanaudenaerde, B. M., and Wuyts, W. A. (2019). Gene correlation network analysis to identify regulatory factors in idiopathic pulmonary fibrosis. Thorax 74, 132–140. doi: 10.1136/thoraxjnl-2018-211929

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Núñez, N. G., Tosello Boari, J., Ramos, R. N., Richer, W., Cagnard, N., Anderfuhren, C. D., et al. (2020). Tumor invasion in draining lymph nodes is associated with Treg accumulation in breast cancer patients. Nat. Commun. 11:3272. doi: 10.1038/s41467-020-17046-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Oshi, M., Asaoka, M., Tokumaru, Y., Angarita, F. A., Yan, L., Matsuyama, R., et al. (2020). Abundance of regulatory T cell (Treg) as a predictive biomarker for neoadjuvant chemotherapy in triple-negative breast cancer. Cancers 12:3038. doi: 10.3390/cancers12103038

PubMed Abstract | CrossRef Full Text | Google Scholar

Pei, Y. Y., Li, G. C., Ran, J., and Wei, F. X. (2017). Kinesin family member 11 contributes to the progression and prognosis of human breast cancer. Oncol. Lett. 14, 6618–6626. doi: 10.3892/ol.2017.7053

PubMed Abstract | CrossRef Full Text | Google Scholar

Picard, E., Verschoor, C. P., Ma, G. W., and Pawelec, G. (2020). Relationships between immune landscapes, genetic subtypes and responses to immunotherapy in colorectal cancer. Front. Immunol. 11:369. doi: 10.3389/fimmu.2020.00369

PubMed Abstract | CrossRef Full Text | Google Scholar

Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson, H., et al. (2019). g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 47, W191–W198. doi: 10.1093/nar/gkz369

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, A. L., Wang, Z. C., De Nicolo, A., Lu, X., Brown, M., Miron, A., et al. (2006). X chromosomal abnormalities in basal-like human breast cancer. Cancer Cell 9, 121–132. doi: 10.1016/j.ccr.2006.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Rojas, L. K., Trilla-Fuertes, L., Gámez-Pozo, A., Chiva, C., Sepúlveda, J., Manso, L., et al. (2019). Proteomics characterisation of central nervous system metastasis biomarkers in triple negative breast cancer. Ecancermedicalscience 13:891. doi: 10.3332/ecancer.2019.891

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, P., Klemp, J. R., Kimler, B. F., Mahnken, J. D., Geier, L. J., Khan, Q. J., et al. (2014). Germline BRCA mutation evaluation in a prospective triple-negative breast cancer registry: implications for hereditary breast and/or ovarian cancer syndrome testing. Breast Cancer Res. Treat. 145, 707–714. doi: 10.1007/s10549-014-2980-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y., Barry, W. T., Seah, D. S., Tung, N. M., Garber, J. E., and Lin, N. U. (2020). Patterns of recurrence and metastasis in BRCA1/BRCA2-associated breast cancers. Cancer 126, 271–280. doi: 10.1002/cncr.32540

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugimoto, Y., Sawant, D. B., Fisk, H. A., Mao, L., Li, C., Chettiar, S., et al. (2017). Novel pyrrolopyrimidines as Mps1/TTK kinase inhibitors for breast cancer. Bioorg. Med. Chem. 25, 2156–2166. doi: 10.1016/j.bmc.2017.02.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Takagi, K., Miki, Y., Shibahara, Y., Nakamura, Y., Ebata, A., Watanabe, M., et al. (2013). BUB1 immunolocalization in breast carcinoma: its nuclear localization as a potent prognostic factor of the patients. Horm. Cancer 4, 92–102. doi: 10.1007/s12672-012-0130-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanaka, H., Phipps, E. A., Wei, T., Wu, X., Goswami, C., Liu, Y., et al. (2018). Altered expression of telomere-associated genes in leukocytes among BRCA1 and BRCA2 carriers. Mol. Carcinog. 57, 567–575. doi: 10.1002/mc.22773

PubMed Abstract | CrossRef Full Text | Google Scholar

Tutt, A., Tovey, H., Cheang, M. C. U., Kernaghan, S., Kilburn, L., Gazinska, P., et al. (2018). Carboplatin in BRCA1/2-mutated and triple-negative breast cancer BRCAness subgroups: the TNT trial. Nat. Med. 24, 628–637. doi: 10.1038/s41591-018-0009-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaclova, T., Gomez-Lopez, G., Setien, F., Bueno, J. M., Macias, J. A., Barroso, A., et al. (2015). DNA repair capacity is impaired in healthy BRCA1 heterozygous mutation carriers. Breast Cancer Res. Treat. 152, 271–282. doi: 10.1007/s10549-015-3459-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidula, N., Dubash, T., Lawrence, M. S., Simoneau, A., Niemierko, A., Blouch, E., et al. (2020). Identification of somatically acquired BRCA1/2 mutations by cfDNA analysis in patients with metastatic breast cancer. Clin. Cancer Res. 26, 4852–4862. doi: 10.1158/1078-0432.ccr-20-0638

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., and Li, M. (2019). Correlate tumor mutation burden with immune signatures in human cancers. BMC Immunol. 20:4. doi: 10.1186/s12865-018-0285-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Weisman, P. S., Ng, C. K., Brogi, E., Eisenberg, R. E., Won, H. H., Piscuoglio, S., et al. (2016). Genetic alterations of triple negative breast cancer by targeted next-generation sequencing and correlation with tumor morphology. Mod. Pathol. 29, 476–488. doi: 10.1038/modpathol.2016.39

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J., Zhang, Y., and Li, M. (2019). Identification of methylation markers and differentially expressed genes with prognostic value in breast cancer. J. Comput. Biol. 26, 1394–1408. doi: 10.1089/cmb.2019.0179

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, C., Gong, J., Jie, Y., Cao, J., Chen, Z., Li, R., et al. (2020). NCAPG is a promising therapeutic target across different tumor types. Front. Pharmacol. 11:387. doi: 10.3389/fphar.2020.00387

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, B., Xu, Y., Woo, J. H., Wang, Y., Bae, Y. K., Yoon, D. S., et al. (2006). Increased expression of mitotic checkpoint genes in breast cancer cells with chromosomal instability. Clin. Cancer Res. 12, 405–410. doi: 10.1158/1078-0432.ccr-05-0903

PubMed Abstract | CrossRef Full Text | Google Scholar

Yusuf, M., Gaskins, J., Mandish, S., May, M. E., Wall, W., Fisher, W., et al. (2020). Tumor infiltrating lymphocyte grade in Merkel cell carcinoma: relationships with clinical factors and independent prognostic value. Acta Oncol. 59, 1409–1415. doi: 10.1080/0284186x.2020.1794033

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Liu, Z., Dong, Y., and Kong, L. (2021). E2F2 drives glioma progression via PI3K/AKT in a PFKFB4-dependent manner. Life Sci. 276:119412. doi: 10.1016/j.lfs.2021.119412

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, E., Li, L., Zhang, W., Wang, W., Chan, Y., You, B., et al. (2019). Comprehensive characterization of immune- and inflammation-associated biomarkers based on multi-omics integration in kidney renal clear cell carcinoma. J. Transl. Med. 17:177. doi: 10.1186/s12967-019-1927-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, S., Ye, Z., and Stanton, R. (2020). Misuse of RPKM or TPM normalization when comparing across samples and sequencing protocols. RNA 26, 903–909. doi: 10.1261/rna.074922.120

PubMed Abstract | CrossRef Full Text | Google Scholar