Calcification-associated molecular traits and therapeutic strategies in hormone receptor-positive HER2-negative breast cancer

Patient samples and study cohorts

We established an HR+/HER2− breast cancer cohort with preoperative mammography (FUSCC-HR+/HER2− mammography). The samples were collected from patients with breast cancer who underwent surgery at Fudan University Shanghai Cancer Center (FUSCC) between January 2013 and December 2014. Detailed clinical information, pathological data, mammography reports, and follow-up data were recorded. Whole-exome sequencing, OncoScan, transcriptomic, and metabolomic data generation, as well as mammography image assessments, were as described in our prior research (IRB ID: 050432-4-1911D)^18,19. Ethical approval for the use of all samples in this study was obtained from the FUSCC Ethics Committee.

Immunohistochemistry (IHC) analysis

Protein expression of specific genes in breast cancer tissues was detected with immunohistochemical staining on paraffin-embedded sections as follows. Deparaffinized slices were immersed in xylene, absolute ethanol, 95% ethanol, and 75% ethanol. The slices were then immersed in sodium citrate antigen retrieval solution and heated over water. Light-blocking was performed with peroxidase blocking solution for 30 min, and was followed by blocking with a 5% BSA solution for 1 h at room temperature. After TBST washing, the diluted antibodies (100 μL) were added dropwise, and the slices were incubated overnight at 4 °C. After TBST washing, corresponding secondary antibodies were added dropwise. DAB working solution was used for color development. Subsequently, the slides were sequentially immersed in 75% ethanol, 95% ethanol, absolute ethanol, and xylene. Finally, neutral resin was used for slide sealing. Immunohistochemistry scoring was determined by positive staining density, which was quantified with a computerized imaging system comprising a Leica charge-coupled device DFC420 camera connected to a Leica DM IRE2 microscope (Leica Microsystems Imaging Solutions). Densities were calculated by counting positive cells within an ×10 high-power field of view (~2 mm²). Representative images of ER, PR, and Ki-67 are shown in Figure S1.

Prediction analysis of microarray 50-gene test classification

The determination of microarray 50-gene test (PAM50) subtypes was based on previously reported methods^20,21. First, mRNA sequencing data were sampled to ensure that the proportion of IHC subtypes in the sample matched that of the training set used for PAM50 subtype determination. Subsequently, the FPKM data were adjusted according to the median expression values of PAM50 genes in the subset described above. Finally, PAM50 classification was determined according to the previously reported criteria.

Estimation of tumor mutational burden

Tumor mutational burden (TMB) was defined as the number of nonsynonymous somatic mutations per megabase (muts/Mb) within the coding region of the captured exome, encompassing missense, nonsense, nonstop, splice site, and translation start site mutations, as well as in-frame and frameshift insertions and deletions.

Estimation of the homologous recombination deficiency score and chromosomal instability score

The homologous recombination deficiency (HRD) score was the sum of the telomeric allelic imbalance (NtAI) score, loss of heterozygosity (LOH) score, and large-scale state transitions (LST) score. The telomeric allelic imbalance score indicated the count of subchromosomal regions exceeding 11 Mb with allelic imbalance extending to the telomere. The LOH score was determined by the number of LOH regions longer than 15 Mb but shorter than the entire chromosome (excluding LOH regions on chromosome 17). The large-scale state transitions score was defined as the count of breakpoints between 2 chromosomal regions longer than 10 Mb after smoothing out regions shorter than 3 Mb. These scores were calculated with copy number data according to previous studies^22–26. The chromosomal instability (CIN) score was derived by summing the squared gene-level GISTIC 2 values²⁷.

Gene set enrichment analysis and single-sample gene set enrichment analysis

Gene set enrichment analysis (GSEA) was conducted to investigate enriched pathways and interpret RNA-seq data, by using predefined gene sets from the Molecular Signatures Database (v.7.1) in GSEA software (v4.0). Default settings were applied for all basic and advanced fields.

Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the characteristic pathway score with the “GSVA” function in the R package “GSVA”²⁸ (v1.42.0). Pathways were downloaded from the Bader laboratory.

Calculation of the sensitivity to endocrine therapy index and CDK4/6 inhibitor sensitivity score

The sensitivity to endocrine therapy index, used to indicate patients’ sensitivity to endocrine therapy, was defined as follows: , where G is the sum of the expression of 18 informative genes (SLC39A6, STC2, CA12, ESR1, PDZK1, NPY1R, CD2, MAPT, QDPR, AZGP1, ABAT, ADCY1, CD3D, CD3D, MRPS30, DNAJC12, DNAJC12, SCUBE2, and KCNE4) and R is the sum of the expression of 10 reference genes (LDHA, ATP5J2, VDAC2, DARS, UGP2, UBE2Z, AK2, WIPF2, APPBP2, and TRIM2)^29–31. Calculation of the CDK4/6 inhibitor sensitivity score was performed with the Gene Set Variation Analysis (GSVA) method with previously reported relevant genes (SAMD2, TIRAP, SETBP1, NNMT, TMSB10, KRT19, THBS2, MDK, EZH2, CDK4, CYR61, SORD, and GSTM3),³² by using the “GSVA” function in the R package “GSVA”²⁸ (v1.42.0).

Metabolic pathway differential abundance score

Differential abundance (DA) scores reflect the tendency for pathways to have elevated levels of metabolites³³. The DA score was calculated as (number of metabolites significantly increased or decreased)/number of measured metabolites within the pathway. DA scores range from −1 to 1. For a given pathway, scores of 1 indicated that all metabolites increased in abundance, whereas scores of −1 indicated that all metabolites decreased in abundance.

Correlation network of metabolites

First, the detected polar metabolites were included for further analysis³⁴. Subsequently, a polar metabolite network was generated by calculation of Spearman’s correlation, and negative correlations were not used. To achieve a concise representation and robust clustering effect, we retained metabolite interactions with correlation > 0.4 and FDR < 0.05, and removed the outlier nodes (KCore < 3) from the metabolite network. The ForceAtlas algorithm was used to lay out the structures for both networks³⁵, and correlation networks were visualized in Gephi (v0.9.3). Specific metabolites are also visualized in Figure 5B. Metabolites in distinct pathways were selected for visualization.

Statistical analysis

The Mann-Whitney Wilcoxon test and Kruskal–Wallis test were applied to analyze continuous variables, whereas Fisher’s exact test was used to compare categorical variables. The P values for comparisons between pairs of groups among the 3 groups were calculated with the Mann-Whitney Wilcoxon test and adjusted for multiple comparisons with the Bonferroni method. Survival curves were constructed with the Kaplan–Meier product limit method and compared with the log-rank test. All analyses were performed in R (v4.1.1) (https://cran.r-project.org/).

FUSCC-HR+/HER2− mammography cohort

To elucidate the molecular features associated with mammographic calcifications in HR+/HER2− breast cancer, we established a cohort of patients with HR+/HER2− breast cancer with preoperative mammography records (Table 1). Primary tumor samples were derived from patients treated at FUSCC between 2013 and 2014. Clinicopathological information and mammography images were collected in detail, in addition to whole exome sequencing (n = 271), OncoScan (n = 199), transcriptomics (n = 312), lipidomics (n = 196), and polar metabolomics (n = 196) data (Figure 1).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Schematic overview of the study design. A multi-omics cohort comprising 316 patients of breast cancer with mammography data. The cohort was stratified according to calcification features. Comparative analyses were conducted across clinicopathological characteristics, multi-omics dimensions, and precision therapeutic strategies.

We assessed the mammography images according to the Breast Imaging Reporting and Data System (BI-RADS)^19,36,37 classification and categorized the calcification status accordingly (Table S1). Calcifications were divided into 3 groups with varying degrees of malignancy suspicion: calcifications assessed as BI-RADS 4C-5 were classified as having high suspicion for malignancy; calcifications assessed as BI-RADS 4A-4B were classified as having low-moderate suspicion for malignancy; and calcifications assessed as BI-RADS 2-3 were classified as probably benign. Tumors without calcifications were categorized as calcification-negative. Due to the limited sample size of patients with calcifications of low-moderate suspicion for malignancy (n = 23), subsequent analyses focused primarily on the remaining patients (Figure 1).

Thus, we established a multi-omics cohort of HR+/HER2-breast cancers with mammography images for further exploration of biological features associated with calcifications, and related implications for precision treatment.

Clinical and pathological characteristics of tumors with different calcification statuses

We first focused on tumors with calcifications of high suspicion for malignancy. These tumors were associated with unfavorable clinicopathological characteristics, and were significantly larger than those in the other groups (in comparison to calcification-negative tumors, P = 0.035; in comparison to tumors with probably benign calcifications, P = 0.061, Mann-Whitney Wilcoxon test) (Figure 2A). These tumors also showed a higher rate of lymph node metastasis and lympho-vascular invasion than tumors without calcifications and tumors with probably benign calcifications (lymph node metastasis rate, 50.5% vs. 44.4% vs. 62.0%, P < 0.001, Fisher’s exact test; lympho-vascular invasion rate, 45.5% vs. 34.9% vs. 58.6%, P = 0.020, Fisher’s exact test) (Figure 2B, C). Moreover, this subset of tumors showed higher Ki-67 IHC staining scores than the other two groups (in comparison to calcification-negative tumors: P = 0.184; in comparison to tumors with probably benign calcifications, P < 0.001, Mann-Whitney Wilcoxon test) (median Ki-67 scores, 20% vs. 20% vs. 30%, P = 0.001, Kruskal-Wallis test) (Figure 2E). Furthermore, these tumors had higher histological grades than the other groups, although the difference did not reach statistical significance (Figure S2).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Clinical characteristics and pathological patterns of different calcification groups. (A) The box plot shows tumor sizes among groups. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 103, 90, and 100, respectively. (B) The stacked bar chart compares the lymph node stages across groups. P values were calculated with Fisher’s exact test. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 103, 90, and 100, respectively. (C) The stacked bar chart shows the lympho-vascular invasion status across groups. P values were calculated with Fisher’s exact test. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 103, 90, and 100, respectively. (D) The box plot illustrates the ages among groups. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 103, 90, and 100, respectively. (E) The box plot shows the immunohistochemistry scores for ER, PR, and Ki-67. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 103, 90, and 100, respectively. (F) Relapse-free survival (RFS) among groups with endocrine therapy. P values were calculated from a two-stage test (probably benign vs. the others). (G) A stacked bar chart illustrates the calcification status for different recurrence outcomes. P values were calculated with Fisher’s exact test.

The patients with probably benign calcifications were significantly older than those with calcifications of high suspicion for malignancy (P = 0.009, Mann-Whitney Wilcoxon test) (Figure 2D), and their tumors exhibited higher ER and PR IHC staining scores than tumors with calcifications of high suspicion for malignancy (P = 0.071 and P < 0.001, respectively, Mann-Whitney Wilcoxon test) (Figure 2E).

Patients with probably benign calcifications tended to have better prognosis with endocrine therapy than the other patients, although this difference did not achieve statistical significance (P = 0.142, two-stage test) (Figure 2F). In addition, there was a higher rate of calcifications with high suspicion for malignancy in patients who relapsed compared to those who did not (46.2% vs. 33.5%, P < 0.001, Fisher’s exact test) (Figure 2G).

Genomic features of tumors among calcification groups

To further understand the molecular features associated with mammographic calcification status, we compared the genomic features among groups. Overall, tumors with calcifications of high suspicion for malignancy had significantly higher tumor mutational burden (TMB), HRD scores, and CIN scores than the other groups (Figure 3A–C). We then analyzed mutational signatures specific to breast cancer and found no statistically significant differences (Figure 3D).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Genomic features of tumors in different calcification groups. (A) The box plot illustrates tumor mutational burden (TMB) scores among calcification groups. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 87, 77, and 87, respectively. (B) The box plot shows homologous recombination deficiency (HRD) scores across groups. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 93, 83, and 88, respectively. (C) The box plot shows the comparison of chromosomal instability (CIN) scores. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 93, 83, and 88, respectively. (D) The proportion of mutation signatures in each sample, on the basis of the Catalog of Somatic Mutations in Cancer (COSMIC) database. (E) Somatic mutations of the top mutated genes. Genes are ordered by total mutation frequency. P values were calculated with Fisher’s exact test (^#P value < 0.1). (F) Somatic CNV among groups, and comparison between tumors with high suspicious calcifications and tumors in the other groups. The P values were calculated with Fisher’s exact test and adjusted for multiple comparisons with the Bonferroni method; the horizontal line in the bottom panel represents the adjusted P value of 0.05. (G) The CNV frequencies of cell cycle genes across groups. P values were determined with Fisher’s exact test.

Mutational profiling indicated that the most frequently mutated cancer-associated genes in the cohort included PIK3CA (46%), TP53 (26%), GATA3 (14%), MAP3K1 (15%), and AKT1 (10%) (Figure 3E). We observed a higher mutation frequency of MAP3K1 in the calcification-negative group than the other groups, but this difference did not reach statistical significance (P = 0.080, Fisher’s exact test) (Table S2).

Furthermore, we systemically characterized somatic copy number variation (CNV) in the FUSCC-HR+/HER2− mammography cohort. Compared to tumors without calcification and those with probably benign calcifications, tumors with calcifications of high suspicion for malignancy displayed deletions in segments of chromosomes 17p (including TP53) and 11q (including ATM) (Figure 3F). Previous reports have suggested that the deletion of chromosome 17q is associated with cell cycle dysregulation³⁸. Targeting the cell cycle is a crucial therapeutic approach in HR+/HER2− breast cancer^4–6,39,40. Hence, we compared the copy number alteration frequency for each cell cycle-associated gene among the groups. Activation of the cell cycle was observed in tumors with calcifications of high suspicion for malignancy (Figure 3G). Notably, genes involved in positive regulation of the cell cycle exhibited an elevated amplification rate in these tumors. For example, the gain or amplification of CCND2 was more pronounced than tumors with probably benign calcifications and those without calcifications (gain/amplification rate, 17.2% in calcification-negative groups vs. 8.4% in probably benign groups vs. 23.9% in high suspicion groups, P = 0.023, Fisher’s exact test). Moreover, genes involved in the negative regulation of the cell cycle showed higher deletion rates in tumors with calcifications of high suspicion for malignancy; the loss or deletion frequency of TP53 was highest among the 3 groups (loss/deletion rate: 39.8% vs. 45.8% vs. 71.6%, P < 0.001, Fisher’s exact test).

In summary, we compared the genomic features of the different calcification groups and observed genomic instability and cell cycle activation in tumors with calcifications of high suspicion for malignancy.

Molecular characteristics and potential therapeutic strategies across distinct calcification groups, as inferred from transcriptomic data

To obtain deeper insights into the transcriptomic characteristics associated with calcifications, we first compared the compositions of established molecular subtypes among tumors differing in calcification status (Figure 4A, B). In the PAM50 subtyping system, the tumors with probably benign calcifications were composed primarily of luminal A and luminal B subtypes, in contrast to tumors without calcifications or with calcifications of high suspicion for malignancy (70.9% in calcification-negative groups vs. 91.0% in probably benign groups vs. 81.6% in high suspicion groups, P = 0.047, Fisher’s exact test) (Figure 4A). In regards to SNF subtype, a taxonomy of HR+/HER2-breast cancer we previously proposed based on this cohort¹⁸, tumors with calcifications of high suspicion for malignancy were mainly SNF3 subtype, in contrast to tumors without calcifications and those with probably benign calcifications (27.1% vs. 29.4% vs. 49.2%, P = 0.064, Fisher’s exact test) (Figure 4B). SNF3 tumors are characterized by activation of the cell cycle and HRD. The higher proportion of the SNF3 subtype was consistent with the higher Ki-67 staining scores (Figure 2E) and HRD scores (Figure 3B) observed in tumors with calcifications of high suspicion for malignancy.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Exploration of transcriptomic features and therapeutic strategies among calcification groups. (A) The stacked bar chart shows the composition of PAM50 subtypes across groups. P values were calculated with Fisher’s exact test. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 103, 89, and 98, respectively. (B) The stacked bar chart compares the composition of SNF subtypes across groups. P values were calculated with Fisher’s exact test. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 59, 51, and 59, respectively. (C) The heatmap shows the comparison of biomarker scores. P values were determined with Mann-Whitney Wilcoxon tests (one group vs. the others) (*P value < 0.05; **P value < 0.01; ***P value < 0.001). (D) The heatmap reveals the estimated abundance of 24 microenvironment cell types. (E) The box plot shows the single-sample gene set enrichment analysis (ssGSEA) scores for ER-associated pathways. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 102, 89, and 96, respectively. (F) The box plots illustrate the sensitivity scores for endocrine therapy among groups. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 102, 89, and 96, respectively. (G) The box plot shows the single-sample gene set enrichment analysis (ssGSEA) scores for cell cycle-associated pathways. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 102, 89, and 96, respectively. (H) The box plots compare the sensitivity scores for CDK4/6 inhibitors among groups. The sample sizes for the calcification-negative, probably benign, and high suspicion groups were 102, 89, and 96, respectively.

Additionally, we curated treatment-associated signatures from the I-SPY2 study to assess molecular characteristics and potential treatment strategies in different groups, on the basis of the RNA data from our cohort (Figure 4C)⁴¹. Each group had unique features. Tumors without calcifications showed a highly immune-activated phenotype, tumors with probably benign calcifications showed a clear ER-associated signature, and tumors with calcifications of high suspicion for malignancy displayed a high degree of proliferation.

Because the tumors among the three groups had distinct transcriptomic characteristics, we next explored the groups separately. Tumors without calcification exhibited immune activation; therefore, we used our RNA data to estimate the relative abundance of 24 microenvironmental cell types among the patients⁴². Tumors without calcification tended to have elevated infiltration of both adaptive and innate immune cells, thus suggesting a hot immune microenvironment (Figure 4D).

For tumors with probably benign calcifications, we performed ssGSEA and confirmed that the upregulation of ER-associated pathways in this subset of tumors, which originated from GSEA^43,44 (Figure 4E, Table S3). A previous study has indicated that ER pathway activation suggests sensitivity to endocrine therapy (ET)⁴⁵. Therefore, we calculated the endocrine therapy index^29–31 for predicting sensitivity to ET. Tumors with probably benign calcifications were found to be more likely to benefit from ET (in comparison to calcification-negative tumors, P = 0.019; in comparison to high suspicion tumors, P < 0.001, Mann-Whitney Wilcoxon test) (Figure 4F).

Tumors with calcifications of high suspicion for malignancy showed higher ssGSEA scores for cell cycle-associated pathways (Figure 4G), thus indicating potential benefits of CDK4/6 inhibitor treatment. We calculated the sensitivity score for CDK4/6 inhibitors³² and found that these tumors were predicted to have elevated sensitivity to CDK4/6 inhibitors (in comparison to calcification-negative tumors, P = 0.007; in comparison to probably benign tumors, P < 0.001, Mann-Whitney Wilcoxon test) (Figure 4H).

In summary, we further explored the molecular characteristics of tumors with various calcification statuses at the transcriptomic level and proposed treatment strategies that may enable calcification-assisted clinical management.

Metabolomic characteristics among calcification groups

Because metabolic reprogramming is a hallmark of cancer⁴⁶, we next explored the metabolic features associated with calcification. We performed metabolic pathway-based DA score analysis to describe the metabolic differences among calcification groups (Figure 5A)⁴⁷. In tumors without calcification, we observed upregulation of vitamin metabolism (particularly thiamine metabolism) and downregulation of lipid metabolism. In tumors with probably benign calcification, thiamine metabolism was among the most down-regulated pathways. We also identified activation of lipid metabolism in these tumors, which has been associated with ER signaling in previous research⁴⁸. In tumors with calcifications of high suspicion for malignancy, upregulation of amino acid metabolism (for example, tryptophan metabolism) was the most prominent metabolic feature. We subsequently conducted network analysis on specific polar metabolites, which indicated similar results (Figure 5B).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Exploratory analysis of metabolomic characteristics among calcification groups. (A) Pathway-based analysis of metabolomic changes among groups. The DA score captures the average gross changes for all metabolites in a pathway. A score of 1 indicates that all measured metabolites in the pathway increase, whereas a score of −1 indicates that all measured metabolites in a pathway decrease. (B) Differentially expressed metabolite correlation network based on polar metabolites with a >0.4 Spearman’s correlation, FDR < 0.05. Correlation networks were partitioned by a graph-clustering algorithm, and the average quantification of groups in the correlation networks is presented.