This study was approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (ID: bc2021063). Patients included in this study were those with breast cancer, with written informed consent, and with successful sequencing results. There are 160 Chinese patients included. Patients were staged according to the TNM system. Samples from these patients were obtained from surgery and fixed in 4% neutral-buffered formalin at 4°C for 24 h. The processed samples were embedded in paraffin wax for storage.

PR, ER, and HER2 expressions were first detected with immunohistochemistry (IHC). Briefly, a 3–5-μm slice of the wax sample was deparaffinized on the glass slide in a 60°C oven for 1 h and washed with xylene for 15 min. The sample was rehydrated with gradient concentrations (100%, 95%, 85%, and 75%) of alcohol. Four minutes of pressure cooking and 10 min of cold-water bath were used to induce antigen retrieval. Hydrogen peroxide (3%) was used to quench endogenous peroxidase activity for 30 min. The slides were blocked with goat serum (10%) for 30 min and stained with rabbit antibodies for ER, PR, and HER2. A secondary antibody biotinylated with horseradish peroxidase (HRP) was further stained for 25 min (K8002 kit, Dako, Santa Clara, CA). The slides were processed with diaminobenzidine for 5 min (K8002 kit, Dako, Santa Clara, CA). After washing with deionized water, hematoxylin was used to restain the slides. ER- and PR-positive rates of >1% were defined as ER+ and PR+, respectively. HER2 score = 3+ and HER2 score = 0/1+ were defined as HER2+ and HER2-, respectively. HER2 score = 2+ was further verified by fluorescence in situ hybridization (FISH). It was performed according to the manual of the PathVysion HER-2 DNA Probe Kit (Abbott Molecular Inc., IL, USA). The HER2 FISH result was interpreted according to the American Society of Clinical Oncology (ASCO)/College of American Pathologists (CAP) guidelines (9).

The tumors were classified into three hormonal subtypes according to ER, PR, and HER2 status. The three hormonal subtypes are the luminal subtype (ER+ or PR+), HER2-enriched subtype (ER- and PR-, HER2+), and TNBC subtype (ER-, PR-, and HER2-).

The library construction and sequencing were performed in the Clinical Laboratory Improvement Amendments (CLIA)/CAP-compliant Molecular Diagnostics Service Laboratory of Shanghai OrigiMed Co., Ltd. Methods were similar to our previous work (10, 11). Briefly, the formalin-fixed paraffin-embedded (FFPE) tissues were deparaffinized by heating in a 60°C oven for 1 h. As a control, white blood cells from the paired whole blood samples were separated by centrifugation.

DNA was extracted using KAPA (Kapa Biosystems, Wilmington, MA, USA) HyperPrep Kit (Illumina, San Diego, CA, USA). Samples with at least 50 ng of double-stranded DNA were applied in further library construction for both tumor tissues and paired white blood cells. A panel of 624 pan-cancer genes was targeted for amplification. To reduce errors from amplification, molecular identifiers (MIDs) were added to DNA segment ends. Barcodes were also added for multiplex sequencing. The prepared DNA libraries were sequenced on an Illumina NovaSeq 6000 sequencer (Illumina, San Diego, CA). Around both ends of the reads, 151 base pairs were read. The average sequencing depth was about 1,000×.

Raw reads were first trimmed for adaptors by Cutadapt (version 1.18). Reads were de-duplicated according to MID labels. The software Burrows–Wheeler alignment with maximal exact matches (BWA-MEM, version 0.7.9a) (12) was used to map all the reads onto University of California, Santa Cruz (UCSC) hg19 reference sequences. The alignment was further recalibrated by BaseRecalibrator of GATK (version 3.8). Short nucleotide mutations were called with Mutect2 (13) and varscan (14). Germline mutations were defined as those in the paired white blood cells but not in the known single-nucleotide polymorphism database (ESP6500, 1000 Genomes, gnomAD, and ExAC). Somatic mutations were retrieved by filtering out possible germline mutations. Germline mutations were obtained from each patient’s white blood sequencing and the single-nucleotide polymorphism databases. Further, copy number variations were called with CNVkit (15), and gene fusions were identified with an in-house pipeline (16).

Tumor mutational burden (TMB) was defined as the number of mutations per million effective coverage length of the genome (1.795717 million bases).

Mutations could be preferentially present in specific genes because of gene length, genomic regions, and other patient characteristics. To reduce mutation bias, MutSigCV (17) was used to identify the significantly mutated genes. The calculation was performed using the service on the GenePattern website (18). Default parameters were used. The resulting P-values were adjusted by the Benjamin and Hochberg method for the correction of multiple testing errors.

A mutation actionable database, OncoKB (19), was used to query the available drugs and their evidence level for each mutation. A python script released by OncoKB, MafAnnotator.py, was used to automate the query. The disease type was set to “BRCA.” According to the evidence level of confidence, actionable mutations were classified into six, namely, “Level_1”, “Level_2A”, “Level_2B”, “Level_3A”, “Level_3B”, and “Level_4”. Level_1 is the treatment with an FDA-recognized biomarker predictive of response to an FDA-approved drug; Level_2A is the treatment with standard-care biomarker predictive of response to an FDA-approved drug. Level_ 2B is the treatment with standard-care biomarker predictive of response to an FDA-approved drug in another indication but not the standard care in this indication. Level_3A is the treatment with compelling clinical evidence that supports the biomarker as being predictive of response to a drug in this indication. Level_3B is the treatment with compelling clinical evidence that supports the biomarker as being predictive of response to a drug in another indication. Level_4 is the treatment with compelling biological evidence that supports the biomarker as being predictive of response to a drug.

A muskoskeletal memorial sloan kettering (MSK) cohort (20) was also used to contrast the results in this study. This cohort included 1,756 patients with breast cancer. Data were downloaded from the cBioPortal (www.cbioportal.org).

Comparisons between categorical variables were performed using Fisher’s exact test. Continuous variables were compared by the Mann–Whitney U test. Multiple-testing comparisons were corrected with the Benjamin and Hochberg method. P-values below 0.05 were considered statistically significant. The similarity score of two factors was obtained by transforming P-values to “-log10(P-value)”. When analyzing the association between clinical characteristics and genomic mutations, the cases were omitted if missing records in the current analysis.

Between 2016 and 2018, 160 patients diagnosed with breast cancer were enrolled in this study ( Table 1 ). The patients were about 50 years old, and 32.5% of them were at late-stage. About 50% of patients were either ER-positive or PR-positive and 28.8% of them were HER2-positive. Samples were harvested from 95 primary tumors (59.4%) and 63 metastatic tumors (39.4%). A family history of cancer (FH) was present in 31.3% of patients. The distribution of FH is shown in Supplementary Figure 1A . Lung cancer, liver cancer, breast cancer, endometrial cancer, and colorectal cancer were diagnosed mostly in family members of patients with breast cancer ( Supplementary Figure 1B ), which accounted for 34.2% of patients with FH.

he median percentage of effective coverage (depth >100) is 99.63% (standard error = 1.56), and the lowest is 90.3%. Short nucleotide genomic mutations in the tumors were called with Mutec2 ( Figure 1 ). Mutations with allele frequency below 0.1% were filtered out for accuracy. The top-10 mutated genes were TP53, PIK3CA, ERBB2, CDK12, PTEN, CCND1, FGF19, FGF3, RB1, and BRCA1 ( Figure 2 ). Only six patients had germline mutations in genes including BRCA1, BRCA2, ATM, and RAD51D. For the top-mutated genes, the association between mutations and clinical characteristics including ER, PR, HER2, and FH is displayed in Figures 2A–D , respectively. CCND1 was found significantly associated with ER status ( Figure 2A ) while ERBB2 and CDK12 were with HER2 status ( Figure 2C ). No significantly associated mutations were found significantly associated with PR ( Figure 2B ) and FH ( Figure 2D ). The association between gene mutations and breast cancer subtypes is shown in Supplementary Figures 2A–C . The HER-enriched subtype had a higher proportion of ERBB2 and CDK12 ( Supplementary Figure 2B ), and the TNBC subtype had a higher proportion of Kirsten Rat Sarcoma Viral Proto-Oncogene (KRAS) and a lower proportion of ERBB2 ( Supplementary Figure 2C ).

The results of this study were contrasted with those of the MSK cohort (20). In the MSK cohort, ER+ and PR+ patients had a higher proportion of PIK3CA, CDH1, MAP3K1, AKT1, ESR1, and GATA3 mutations and a lower proportion of TP53 mutations ( Supplementary Figures 2D, E ); HER2+ patients had a higher proportion of TP53, NF1, and MLL2 and a lower proportion of CDH1, GATA3, and MAP3K1 ( Supplementary Figure 2F ). The luminal subtype had a higher proportion of PIK3CA, CDH1, GATA3, and MAP3K1 mutations and a lower proportion of TP53 mutations ( Supplementary Figure 2G ), and the HER2-enriched subtype had a higher proportion of TP53 and NF1 mutations and a lower proportion of GATA3 mutations ( Supplementary Figure 2H ). No significant mutations were found in the TNBC subtype ( Supplementary Figure 2I ).

To reduce mutational bias from gene length, genomic region, and patient characteristics, etc., MutSigCV was used to identify the significantly mutated genes. Five genes (TP53, PIK3CA, AKT1, PTEN, and GATA3) were found significantly mutated with an adjusted P-value less than 0.05. None of these genes was significantly associated with ER, PR, HER2, and FH.

Testing each gene's mutation frequency by Fisher's exact test had a low power to find out significant genes because of multiple testing problem. To improve statistical power, a clustering method was employed in advance. The clustering procedure included three steps. First, the association between two mutations was calculated with Fisher’s exact test. Second, the resulting P-values were transformed into a similarity score between two genes as “-log10(P-value)”. Third, the similarity scores were used to hierarchically cluster genes with Euclidean distance and complete linkage ( Figure 3 ). From the heatmap result, four clusters having a high association within them were marked with rectangles. Patients were classified into four corresponding groups if there is a gene mutation in each gene cluster. Mutations within each group of patients were annotated with KEGG pathways. A hypergeometrical test was used to infer the significance of enrichment ( Figure 4 ). Interestingly, group 1 was enriched with signaling pathways regulating pluripotency of stem cells (21) ( Figure 4A ); group 2 was enriched with protein digestion and absorption pathway and extracellular matrix (ECM)–receptor interaction pathway (22) ( Figure 4B ); group 3 was enriched with hippo signaling pathway (23) ( Figure 4C ); and group 4 was enriched with VEGF signaling pathway (24) ( Figure 4D ).

The transformed P-value was calculated after Fisher’s exact test between the status of clinical characteristics and each group of patients ( Figure 5A ). Group 1 was significantly associated with PR status; groups 2 and 3 were significantly associated with ER status. Fisher’s exact test was also used to compare the gene cluster and FH status, which indicated that cluster 2 and cluster 4 were found negatively associated with a family history of cancer and age, respectively ( Figure 5B ). The association of breast cancer subtypes and TNM staging with mutation clusters is shown in Supplementary Figures 3A, B . No significant associations were found between them. We also analyzed the association between gene mutations and clinical characteristics in the MSK cohort. Cluster 2 showed a significant association with PR status.

Mutations in the tumors of the four groups of patients were annotated with the OncoKB database (19). Actionability was classified into six levels according to the confidence in the evidence ( Figure 6 ). Level_1 was of the highest evidence, and Level_4 was of the lowest evidence. The highest evidence level, Level_1, was the treatments with an FDA-recognized biomarker predictive of response to an FDA-approved drug, and the lowest evidence level, Level_4, was the treatments with compelling biological evidence that supports the biomarker as being predictive of response to a drug. Patient groups 1, 2, 3, and 4 had 24.1%, 36.5%, 38.7%, and 41.3% of patients with an FDA-recognized biomarker predictive of response to an FDA-approved drug, respectively ( Figure 6 ). ERBB2 and PIK3CA were the most enriched actionable genes in groups 2, 3, and 4 with a high evidence level. Although the percentage of actionable mutations in FGFR1 was high in group 1, they had a low evidence level for majority of the patients.