Perturbation and stability of PAM50 subtyping in population-based primary invasive breast cancer

Abstract

PAM50 gene expression subtypes represent a cornerstone in the molecular classification of breast cancer and are included in risk prediction models to guide therapy. We aimed to illustrate the impact of included genes and biological processes on subtyping while considering a tumor's underlying clinical subgroup defined by ER, PR, and HER2 status. To do this we used a population-representative and clinically annotated early-stage breast tumor cohort of 6233 samples profiled by RNA sequencing and applied a perturbation strategy of excluding co-expressed genes (gene sets). We demonstrate how PAM50 nearest-centroid classification depends on biological processes present across, but also within, ER/PR/HER2 subgroups and PAM50 subtypes themselves. Our analysis highlights several key aspects of PAM50 classification. Firstly, we demonstrate the tight connection between a tumor's nearest and second-nearest PAM50 centroid. Additionally, we show that the second-best subtype is associated with overall survival in ER-positive, HER2-negative, and node-negative disease. We also note that ERBB2 expression has little impact on PAM50 classification in HER2-positive disease regardless of ER status and that the Basal subtype is highly stable in contrast to the Normal subtype. Improved consciousness of the commonly used PAM50 subtyping scheme will aid in our understanding and interpretation of breast tumors that have seemingly conflicting PAM50 classification when compared to clinical biomarkers. Finally, our study adds further support in challenging the common misconception that PAM50 subtypes are distinct classes by illustrating that PAM50 subtypes in tumors represent a continuum with prognostic implications.

Conflict of interest statement

The authors declare no competing interests.

© 2023. Springer Nature Limited.

Figures

Fig. 1. Patterns of PAM50NC versus PAM50NC_2nd subtype.

In panels (a–e), the left panels show the cross-tabulated PAM50NC subtype versus the PAM50NC_2nd subtype for separate tumor subsets, whereas the right panels show the corresponding difference (delta) in Spearman correlation between PAM50NC and PAM50NC_2nd subtype based on the average Spearman correlation of the 100 NC classifications for each case. In the cross tables, colored boxes highlight consistent subtype patterns between PAM50NC and PAM50NC_2nd subtypes. Of all 6233 tumors, 6228 had an unambiguous second-best subtype based on NC classification. a All SCAN-B tumors. b TNBC tumors. c ERnHER2p tumors. d ERpHER2p tumors. e ERpHER2n tumors. f Heatmap of Pearson correlations between PAM50 centroids. Heatmap cells marked with colored boxes show centroid correlation patterns consistent with the PAM50NC and PAM50NC_2nd subtype patterns shown in panels (a–e). g Scatter plot of LumA correlation values versus LumB correlation values for tumors classified as LumANC – LumBNC_2nd or LumBNC – LumANC_2nd (n = 1599). The red line corresponds to a 1:1 relationship between correlation estimates. Boxplot elements correspond to: (1) center line = median, (2) box limits = upper and lower quartiles, (3) whiskers = 1.5x interquartile range.

Fig. 2. Association of PAM50NC_2nd subtype with patient outcome.

a Overall survival (OS) for endocrine-treated ERpHER2nLNn patients >50 years of age that were PAM50NC subtyped as LumANC. Patients are stratified by their PAM50NC_2nd subtype. b Overall survival for endocrine-treated ERpHER2nLNn patients >50 years of age that were PAM50NC subtyped as LumBNC. Patients are stratified by their PAM50NC_2nd subtype. c Distributions for age at diagnosis (left), tumor size (center), and ROR T0 scores (right) obtained from ref. in endocrine-treated ERpHER2nLNn patients >50 years of age comparing cases subtyped as LumANC – LumBNC_2nd versus LumANC – NormalNC_2nd. d Distributions for age at diagnosis (left), tumor size (center), and ROR T0 scores (right) obtained from ref. in endocrine-treated ERpHER2nLNn patients >50 years of age comparing cases subtyped as LumBNC – HER2ENC_2nd versus LumBNC – LumANC_2nd. e Overall survival for endocrine-treated ERpHER2nLNn patients >50 years of age that were PAM50NC subtyped as LumANC and as ROR-low risk category according to ref. . Patients are stratified by their PAM50NC_2nd subtype. f Overall survival for endocrine-treated ERpHER2nLNn patients >50 years of age that were PAM50NC subtyped as LumBNC and as ROR-high risk category according to ref. . Patients are stratified by their PAM50NC_2nd subtype. Boxplot elements correspond to: (1) center line = median, (2) box limits = upper and lower quartiles, (3) whiskers = 1.5x interquartile range.

Fig. 3. Study overview and PAM50 reclassification results for the leave-oneGeneCluster-out strategy.

a Study outline, perturbation methodology, and subtype switch concept. A sample is called as having a subtype switch if the PAM50NC subtype is observed in ≤50% of the 100 PAM50perturb reclassifications (right panel). b Left panel, size of identified SRIQ core gene clusters defined from 9206 RNA sequencing profiles from ref. . Center panel, heatmap of average PAM50 centroid value for each gene set for each PAM50 centroid subtype. Right panel, Spearman correlation of average SRIQ FPKM gene cluster expression for each gene set combination in all 9206 RNA sequencing profiles. c Spearman correlation matrix of average SRIQ FPKM gene cluster expression versus rank-based scores for eight reported biological metagenes from Fredlund et al. for the 6233 tumors included in this study. d Heatmap of scaled FPKM expression for PAM50 genes stratified by SRIQ gene cluster definition and ordered by clinical group and PAM50NC subtype for the 6233 included tumors. e Percent of tumors switching subtype (i.e., a different PAM50perturb subtype compared to PAM50NC) by the leave-oneGeneCluster-out strategy on a whole cohort level stratified by PAM50NC subtypes for the 6233 included tumors. f Percent of tumors switching subtype by the leave-oneGeneCluster-out strategy on a whole cohort level stratified by tumors’ ER, PR, and HER2 status.

Fig. 4. PAM50 reclassification results for the leave-oneGeneCluster-out strategy when stratified for molecular and clinical subgroup.

a Heatmap showing the proportion of tumors that switched subtype after gene set exclusion stratified by molecular and clinical subgroup. Numbers represent total group sizes per row. b Top panel shows the proportion of ERpHER2n LumANC tumors with a PAM50perturb subtype different from their PAM50NC subtype, i.e., switching subtype, when excluding a specific gene set in leave-oneGeneCluster-out reclassification. Lower panel shows the distribution of the PAM50perturb subtypes in tumors that switched in the top panel, numbers on top represent the total number of samples that switched subtype. c The same illustration as in (b), but for ERpHER2n LumBNC tumors. d Heatmap showing the proportion of tumors that switched subtype after gene set exclusion that had a PAM50perturb subtype similar to the PAM50NC_2nd subtype. Numbers represent total group sizes per row. e Summary bar plots of the percentage of tumors in each clinical group further stratified by their PAM50NC subtype that never switch subtype across all gene set perturbations, i.e., the PAM50perturb subtype is the same as the PAM50NC subtype in all perturbations. These cases are hereon referred to as K0 cases.

Fig. 5. Leave-oneGeneCluster-out perturbation and association to patient outcome.

Forest plot of hazard ratios with 95% confidence intervals from univariate Cox regression, using DRFI as clinical endpoint, for tumors that switched subtype versus tumors that did not switch subtype (reference) after exclusion of a gene set in a TNBC tumors, b ERpHER2n tumors, and c endocrine-treated ERpHER2n tumors only. d Kaplan–Meier plot of DRFI for PAM50perturb subtypes in endocrine-treated ERpHER2n LumANC tumors after exclusion of gene set 1 (proliferation). e Kaplan–Meier plot of DRFI for PAM50perturb subtypes in endocrine-treated ERpHER2n LumBNC tumors after exclusion of gene set 3 (basal keratins). f Boxplots of rank-based scores for the mitotic progression, basal, steroid response, and lipid metagenes for endocrine-treated ERpHER2n LumANC tumors in panel (d). g Boxplots of rank-based scores for the mitotic progression, basal, steroid response, and lipid metagenes for endocrine-treated ERpHER2n LumBNC tumors in panel (e). Note that not all included cases in the study have DRFI outcome data, thus the difference in sample numbers between boxplots and survival plots. Boxplot elements correspond to: (1) center line = median, (2) box limits = upper and lower quartiles, (3) whiskers = 1.5x interquartile range.

Fig. 6. Refined single sample PAM50 subtyping in ERpHER2n tumors based on leave-oneGeneCluster-out perturbation stable tumors.

a Outline of the scheme to create refined ERpHER2n PAM50 centroids (termed PAM50K0) used for single sample classification by Spearman correlation based on FPKM values only (i.e., no gene centering). b Sankey plot of subtype change for ERpHER2n tumors when performing PAM50K0 classification as outlined in (a). c Kaplan–Meier plot of DRFI for PAM50K0 subtypes in endocrine-treated ERpHER2n tumors. d Boxplots of rank-based scores for the mitotic checkpoint, steroid response, and basal metagenes for endocrine-treated ERpHER2n tumors stratified by PAM50K0 subtypes. e Kaplan–Meier plot of DRFI for PAM50K0 subtypes in endocrine-treated ERpHER2n LumANC tumors. HER2EK0 and NormalK0 groups excluded due to size. f Left panel, Kaplan–Meier plot of OS for PAM50K0 subtypes in all endocrine-treated ERpHER2n LumANC tumors. HER2EK0 and NormalK0 groups excluded due to size. Right panel, same plot but only for non-K0 tumors (i.e., tumors not included in the PAM50K0 centroid creation). g Boxplots of rank-based scores for the mitotic checkpoint, steroid response, and basal metagenes for endocrine-treated ERpHER2n LumANC tumors stratified by PAM50K0 subtypes. Note that not all included cases in the study have DRFI outcome data, thus the difference in sample numbers between boxplots and survival plots. Boxplot elements correspond to: (1) center line = median, (2) box limits = upper and lower quartiles, (3) whiskers = 1.5x interquartile range.