Genetic landscape of breast cancer and mutation tracking with circulating tumor DNA in Chinese women

Patients’ characteristics

A total of 273 breast cancer patients were enrolled in the study (Supplementary Table 1). All patients were women, and the median age was 52 years. The patients were categorized based on tumor subtype, with 77 luminal A (28.2%), 74 luminal B (27.1%), 67 HER2+ (24.5%), and 45 basal-like (16.5%) patients included. Some of the patients had received monotherapy with an aromatase inhibitor (AI), while others had been administered combination AI or fulvestrant and CDK4/6 inhibitor therapy. However, no patients had received prior poly ADP ribose polymerase (PARP) inhibitor therapy.

The ctDNA mutational landscape in Chinese women with breast cancer

High-frequency molecular variants were analyzed in pairs of tumor tissue and plasma. The consistency, defined as harboring at least one common gene mutation site, was evaluated between the tumor tissue and paired plasma samples. The consistency values between tissue and plasma samples in stage I, II, III, and IV breast cancers were 80% (32/40), 80% (36/45), 100% (13/13), and 100% (3/3), respectively. The top three most prevalent mutated genes were as follows: PIK3CA (tissue: 42.31%, plasma: 24.04%), TP53 (tissue: 32.69%, plasma: 43.27%), and ERBB2 (tissue: 19.23%, plasma: 14.42%). The top three variants were PIK3CA p.H1047R (c.3140A>G), PIK3CA p.E545K (c.1633G>A), and PIK3CA p.E542K (c.1624G>A) (Figure 1). Our findings were consistent with those reported by Andre et al [20].

Next, matched DMI-tagged cell-free DNA (cfDNA) was sequenced using the Omigen 101 genes panel kit at a targeted sequencing depth of 40,000×–50,000×. The mutation profile was derived from 205 detected female Chinese patients’ cfDNA with an overall mutation detection rate of 94.15% (193/205; Figure 2). The allele fraction (AF) detected in our cohort ranged from 0.05%–32.48%. Some of these mutations were actionable, including SNVs, such as PIK3CA p.H1047L, p.E545K, and p.E542K, BRCA1/2 loss-of-function oncogenic mutation, and copy number mutations, such as ERBB2 amplification. No mutation was detected from the panel of the remaining 12 patients. Moreover, the patients with negative cfDNA tests were primarily in stage I or II of disease, including six patients in stage I and four in stage II (Figure 2).

These data suggest that somatic mutations with frequencies as low as 0.05% could be detected in DMI-tagged cfDNA of early-stage breast cancer patients with a targeted sequencing depth of 40,000×–50,000×.

Mutation landscape across different molecular subtypes

The distribution of mutations in patients was further analyzed based on four molecular subtypes. The prevalence of changes in cancer-related genes/pathways was compatible with that reported in previous studies of treated advanced breast malignancies, including frequent oncogenic mutations in the PI3K pathway and loss-of-function mutations in the DNA damage response and tumor suppressor pathways [21, 22]. Frequently varied genes in all four molecular subtypes are presented in Figure 3. The PIK3CA detection rates in basal-like (38.24%, 13/34) were higher than in luminal B (31.67, 19/60), HER2+ (30.61%, 15/49), and luminal A (26.98 %, 17/63) breast cancers (Figure 3A). Meanwhile, PTEN mutations were more frequently detected in HER2+ and basal-like subtypes and were less common in the luminal type, with detection rates of 14.29% (7/49), 11.76% (4/34), 9.52% (6/63), and 8.33% (5/60), respectively (Figure 3A). BRCA1/2 mutations were predominantly detected in basal-like patients (14.71%, 11.76%; Figure 3B), while TP53 mutations were frequently detected in all four molecular subtypes (38.24%, 36.73%, 31.75%, and 25% for basal-like, HER2+, luminal A, and luminal B, respectively; Figure 3C).

cfDNA yield and ctF (highest cfDNA allele frequency) correlate with cancer stage

To investigate the correlations between cfDNA and corresponding clinical characteristics, we analyzed cfDNA yield in healthy participants and breast cancer patients with different disease stages. Considering that tumors generate a certain amount of free nucleotides into the peripheral blood, the cfDNA yield could represent the amount of released tumor DNA in the patient’s body [23, 24]. We found that cfDNA was significantly associated with disease stage (p = 0.033, r = 0.9; Figure 4A). Meanwhile, ctF, defined as the highest cfDNA allele frequency in patients with multiple alterations, was found to be significantly associated with disease stage (p = 0.032, r = 0.97; Figure 4B).

Number of ctDNA mutations in different molecular subtypes and disease stages

Cancers with higher tumor mutational burdens (TMBs) are at a higher risk of being recognized by immune cells [25]. Hence, we used the total number of ctDNA non-synonymous coding mutations detected per megabase to calculate TMB and detected the highest TMB in the basal-like subtype, achieving statistical significance in comparison with the luminal subtypes (luminal A: p = 0.0047, luminal B: p = 0.044; Figure 5B). Our data also shows that the total number of ctDNA mutations positively correlates with the breast cancer stage (p = 0.114, r = 0.89; Figure 5A), although this result was not statistically significant.

NGS detection using tissue samples better reflects ERBB2 copy number amplification than using blood samples

We next sought to explore the application of NGS detection for predicting ERBB2 copy number amplification in breast cancer patients. To this end, we applied a 101 gene panel for detecting ERBB2 copy number amplification in 123 tissues and 202 blood samples. Immunohistochemistry (IHC) and Fluorescence in situ hybridization (FISH) were performed simultaneously. We then analyzed and compared the ERBB2 copy number amplification results obtained by all the methods. The NGS results for tissue samples were highly consistent with IHC and FISH, with a sensitivity of 86%, a specificity of 99%, and accuracy of 95% (ratio of ERBB2 duplication region ≥75%, The duplication of ERBB2 is identified when the proportion of amplified ERBB2's exon region greater than or equal to 0.75). However, the results of NGS for blood samples were relatively inconsistent with those of IHC and FISH results, presenting a sensitivity of only 5.4%, a specificity of 97%, and accuracy of 71% (Supplementary Table 2).

For samples with an IHC test result of 2+, an additional FISH test was required to determine the final ERBB2 expression. We analyzed the FISH and NGS results for tissue samples simultaneously to evaluate the consistency between the two detection methods in IHC 2+ samples. We obtained 100% consistency in 12 tissue samples, implying the perfect accuracy of NGS for the detection of ERBB2 copy number amplification in IHC 2+ samples (Supplementary Table 3). Our findings, therefore, support the application of NGS for the detection of ERBB2 copy number amplification in tissue samples, particularly when combined with NGS and IHC, which would serve to improve the accuracy of detecting ERBB2 copy number amplification.

Dynamic monitoring and prognosis assessment

The application of ctDNA for dynamic monitoring and prognostic assessment of tumor patients was also explored. We focused on serial samples of 15 patients who underwent two or more ctDNA tests.

Eleven patients presented with decreased ctF ratio (ctF ratio: ctF2/ctF1, ctF1 refers to the highest AF in the first test, and ctF2 refers to the highest AF in the second test) (Figure 6, Supplementary Table 4), which is consistent with their imaging findings: 10 patients had significantly shrunken tumor size, and 1 patient had no significant increase in tumor size. Meanwhile, four patients had increased ctF ratios during or after treatment, 3 of whom P112, P219, and P86, had metastasis or relapse. The remaining patient required continued follow-up and observation to determine the clinical outcome.

Patient P112 had locally advanced breast cancer, which was first diagnosed in December 2015, and received chemotherapy, endocrine therapy, and radiotherapy. The disease then progressed in March 2017 with liver metastasis. The first ctDNA detection (T1) was performed in April 2018, and postoperative lumbar spine metastasis occurred during treatment. The second detection (T2) was performed four months later when the disease had continuously progressed. Based on the patient’s genetic detection results, we recommended alpelisib plus fulvestrant treatment. We tracked the PIK3CA p.E545K and TP53 p.E286K mutations and found that their frequency increased from 0.64% and 0.18% to 2.65% and 1.54%, respectively (Figure 7). The change in ctDNA in the patient’s plasma was consistent with the clinical disease progression, resulting in a poor prognosis. A similar situation occurred in patient P219 with late-stage breast cancer, and the TP53 p.Y220C mutation increased from 0.16% to 0.97%, with a ctF ratio of 2.310. Additionally, patient p311 was diagnosed with HER2+ breast cancer and received neoadjuvant chemotherapy and target therapy. The patient presented with decreased ctF ratio and tumor size during mutation tracking and reached complete clinical response during surgery.

To further understand the molecular characteristics and clonal evolution in our study samples, clonal analysis was performed for several patients using mutation tracking data. Certain clones were undetectable following treatment in patients P080, P105, P268, and P391. Moreover, ctF values decreased after receiving treatment in these patients, suggesting a decreased tumor mutational burden or disappearance of specific clonal clusters resulting in a relatively better prognosis (Figure 8). However, new emerging clones were observed following treatment in patient P086.

Patient P086 was diagnosed with triple-negative breast cancer (TNBC) in April 2018 and received the first ctDNA detection (T1) before standard neoadjuvant chemotherapy with no significant change observed in the tumor during this period. The patient then discontinued our treatment schedule for four months and came back with a larger mass. The second ctDNA detection was performed in February 2019 (T2), after which an additional four rounds of chemotherapy were performed. No significant change in the tumor was observed during the treatment period. The third ctDNA detection and surgery were performed in May 2019. The newly emerged mutation at T2 was ERBB2 p.R678Q, which was also detected at T3 with an increase in mutation frequency from 2.81% to 3.9%. The ctF value then increased, while the clone analysis demonstrated an increase in clones. These analyses indicated that the patient exhibited a tendency for disease progression, consistent with the patient’s clinical manifestations. Studies have shown that liquid biopsy is more informative than traditional methods, such as the CA-125 test and imaging, capable of predicting tumor recurrence 7–10 months earlier [26].

Moreover, previous reports indicate that ERBB2 p.R678Q likely represents an oncogenic mutation [27]. Herein, the cytological study results showed that anti-HER antibody drugs, trastuzumab and pertuzumab, as well as small molecule kinase inhibitors, lenatinib, and lapatinib, elicited inhibitory effects on tumor cells harboring ERBB2 p.R678Q mutations [28]. Although patient P086 tested negative for HER2 amplification using IHC and FISH, an oncogenic missense mutation in the HER2 subtype was detected using ctDNA, suggesting that HER2-targeted therapy may be an effective follow-up treatment option.

Considering that patient P086 was relatively insensitive to standard neoadjuvant chemotherapy regimen, we next sought to analyze the HRD (homologous recombination deficiency) status of several TNBC patients, including patient P086, to explore optional individualized treatments for TNBC patients with different HRD statuses. The HRD score predicts TNBC or BRCA1/2-mutated patient response to DNA damaging neoadjuvant chemotherapy. Patients with an HRD score ≥ 42 were defined as HR deficient and tended to be more sensitive to platinum-containing regimens [29]. Meanwhile, sequencing data analysis revealed that the HRD-sum of patient P086 was remarkably high (Supplementary Table 5), which might be related to the patient’s relatively poor prognosis. Among the five TNBC patients, two were HR deficient, accounting for 40% of the tested patients; the result was consistent with that reported in previous literature (42%) [30]. HRD-positive patients might not benefit from conventional chemotherapy, neoadjuvant chemotherapy, or targeted therapy; however, PARP inhibitors may represent an effective treatment option for these patients [31].

Patients and samples

Thirty healthy and 273 breast cancer patients were enrolled at The First Affiliated Hospital of Wenzhou Medical University. All patient-derived specimens and patients’ information were collected and archived under protocols approved by the institutional review board of the First Affiliated Hospital of Wenzhou Medical University. Informed consent was obtained from each patient for the use of their blood and tissue samples. All patients underwent complete tumor staging according to the seventh edition tumor, node, and metastasis (TNM) criteria of breast cancer [36]. A 10-mL sample of whole peripheral blood was collected at initial diagnosis in cfDNA blood collection tubes (Omigen, Hangzhou, China) and was processed for plasma and white blood cells (WBC) within 72 h, resulting in the collection of approximately 4 mL of plasma. The first blood sample from each patient was collected before tissue biopsy/surgery, and some patients had multiple time points of follow-up sampling after therapy. Tumor tissues were collected and either analyzed fresh or were formalin-fixed and paraffin-embedded (FFPE). Her2 status was tested using IHC. Samples assigned a value of 2+ were further confirmed using FISH. From April 2018 to September 2019, there were 205 cfDNA and 131 tissue/FFPE samples; 104 samples were paired with tissue/FFPE and cfDNA samples.

DNA extraction

Genomic DNA was extracted from the tumor tissue, FFPE, or WBC using an Omigen extraction kit (Omigen, Inc., Hangzhou, China), followed by fragmentation using an Omigen fragmentation kit (Omigen, Inc., Hangzhou, China). cfDNA was extracted from plasma using an OMInano isolation kit (Omigen, Inc., Hangzhou, China). DNA yield was quantified with Qubit HS DNA (Life Technologies, Carlsbad, CA, USA). All the procedures were performed according to the manufacturer’s instructions.

Library construction

A DNA library was prepared using an Omigen lib prep kit (Omigen, Inc., Hangzhou, China). Briefly, DNA was end prepared, then ligated with barcoded molecular adaptors, followed by polymerase chain reaction (PCR) amplification to add indexes. Library yield was quantified with Qubit, and the size was measured using a 2100 bioanalyzer (Agilent, Santa Clara, CA, USA).

Hybrid capture

To ensure minimal coverage, libraries for each pool were prepared according to the Omigen libraries pooling guideline (China patent: 201910515833.4). Libraries were pooled and target enriched using the Omigen 101 genes panel (Supplementary Table 6) and Omigen hybrid capture kit (Omigen, Inc., Hangzhou, China). The captured pool was then amplified using PCR, and quality control (QC) was performed using a 2100 bioanalyzer. The enriched pool was sequenced using Illumina X10 (Illumina, Inc., USA) with paired-end reads and average sequencing depth of 8,000×–13,000× and 40,000×–50,000× for tissue DNA and DMI-tagged cfDNA samples, respectively.

Sequence data processing and analysis

Somatic single-nucleotide variant (SNV), fusion, and copy number variation (CNV) were analyzed and annotated using the Omigen 101 genes panel (Supplementary Table 6) and Omigen bioinformatics program pipeline including QC, alignment, and variant annotation, followed by filtering. Thirty cfDNA samples from healthy participants were used as background controls.

Exome FASTQ sequences were mapped to the human assembly NCBI build 37 (hg19) using the Bowtie2 alignment software [37]. SAMtools (v0.1.19) was used to convert SAM files to compressed BAM files and sort the BAM files using chromosomal coordinates [38]. PCR duplicates were marked using Picard (v1.131) and subsequently sorted. Marked BAM files were realigned using a Genome Analysis Toolkit (GATK-3.6) at intervals with indel mismatches. GATK was used for performing local realignment for indel calls, and VarDict was used for detecting SNVs [39, 40]. Germline variants were removed by filtering out variants in normal blood samples. Gene annotations and function prediction scores were performed using Annovar, dbSNP build135, 1000Genomes, Polyphen, Avsift, and COSMIC [41–44]. Non-synonymous somatic mutations were identified in the bulk exome sequencing data of five TNBC patients with two or three longitudinal samples. Bulk tumor genomic copy number profiles were estimated from the pair-end exome sequencing depth using the CNV kit [45]. Mutation clusters with only one mutation were excluded from further analysis. The PyClone cluster frequencies were calculated as the mean variant allele frequencies of mutations within each cluster and were visualized using the R package ‘clonevol’ [46, 47].

HRD score calculation

A portion of sporadic TNBCs and BRCA1/2-mutated tumors have DNA repair defects and are sensitive to DNA damage treatment [31]. Based on the loss of heterozygosity (LOH), telomere allelic imbalance (TAI), and large-scale state transition (LST), three independent DNA-based genomic instability measurements were developed. The HRD score represents the arithmetic sum of LOH (number of LOH regions >15 Mb but less than the length of a whole chromosome) + TAI (regions of allelic imbalance that extend to the subtelomere but do not cross the centromere) + LST (breakpoints between regions of imbalance >10 Mb after filtering out regions <3 Mb). HRD scores range between 0 and 100. Tumors with HRD scores ≥42 or BRCA1/2 mutations were defined as having an HR deficient status. Tumors with an HRD score of <42 that also lacked the BRCA1/2 mutation were defined as HR non-deficient. The threshold of 42 was selected based on the fifth percentile of HRD scores in tumors with known BRCA1/2 mutations or methylation status [48, 49]. This process was achieved using a scarHRD software based on whole-exome sequencing data [50].

Statistical analyses

A t-test was performed to compare ctDNA yield, mutation number, and the highest AF between different disease stages and different molecular subtypes. Moreover, their correlation was tested using Pearson’s correlation coefficient. All P values were two-sided, and values <0.05 were considered statistically significant. All statistical analyses were performed using the R language, version 3.6.0.