Proteogenomic analysis dissects early-onset breast cancer patients with prognostic relevance
Breast cancer is the most common malignancy among women globally and is the leading cancer among women in Korea . The occurrence of young breast cancer (YBC), i.e., breast cancer in patients aged 40 years or younger at diagnosis, is rare; however, its proportional incidence in Korea is approximately 10% , which is twice the incidence rate in the US . Patients with YBC have a genetic susceptibility caused by BRCA1/2 pathogenic variants, showing aggressive clinical phenotypes with advanced stages, high-grade and triple-negative breast cancer (TNBC), and an unfavorable prognosis, including late recurrence – . However, the molecular features of YBC underlying these more aggressive tumor characteristics and worse clinical outcomes are unclear. Several studies have revealed the characteristics of YBC. With respect to somatic mutations, the Cancer Genome Atlas (TCGA) and Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) databases revealed that premenopausal (PreM) patients with breast cancer have greater numbers of mutations in 5 genes— CDH1 , GATA3 , MLL3 , GPS2 , and PI3KCA . Moreover, compared to post-menopausal (PostM) tumors, gene expression in PreM tumors is enriched in the integrin and laminin signaling pathways, EGFR signaling activation, and TGF-β, especially in estrogen receptor-positive (ER + ) breast cancer . Another study revealed that, compared with Western patients, Asian PreM patients presented a characteristic increase in tumor-infiltrating lymphocytes (TILs) and a decrease in TGF-β signaling, which suggested that Asian YBCs harbor a more immune-active microenvironment . An age-specific gene expression study of breast cancer in Middle Eastern women identified 63 genes that were specific to YBC ( ≤ 45 years) and 2 genes ( TIAM1 and VANGL2 ) whose expression was significantly lower in breast tumors of very young ( ≤ 35 years) women than in breast tumors from other age groups . In recent studies, the application of proteogenomic (PG) analysis has improved our understanding of breast cancer by providing comprehensive molecular signatures relevant to clinical features – . Stromal-enriched clusters, G protein-coupled receptor clusters, and several targetable biological pathways were also identified. Although this approach revealed new clusters and refined subgroups with potential clinical benefits in breast cancer, there is a need to investigate the molecular landscape in young Asian patients. Therefore, this study aimed to evaluate and analyze the molecular characteristics of YBC in Korean patients (≤40 years) via integrated multiomics studies, including genomic, transcriptomic, and proteomic features, and to identify pathways for predicting the clinical outcomes of patients with YBC. 

Study objectives and specimens A total of 178 participants with histologically defined breast cancer, aged 40 years or younger, treated at the National Cancer Center (NCC) in Korea were included in this study. Tumor and adjacent normal tissue samples were obtained from surgically resected specimens, and blood samples were collected from the patients. This study was approved by the Institutional Review Board of the NCC (IRB nos. NCCNCS 13717, NCC2017-0127, and NCC2020-0135). Since some patients had limited samples, only 126 patients for whom all proteogenomic data were fully generated were used for analysis. We retrospectively reviewed the medical and pathological records of the patients, including histological diagnoses of surgical specimens, tumor staging, histological grade, treatment history (type of surgery, use of chemotherapy, hormone therapy, anti-HER2 therapy, and radiotherapy), and follow-up data (recurrence, metastasis, and death). Tissue specimens and blood samples were provided by the NCC Bio Bank of the National Cancer Center, Korea. The frozen tissue samples were weighed and washed with cold phosphate-buffered saline (PBS) to remove blood contamination, placed in a tube (Covaris, Woburn, MA, USA), snap-frozen in liquid nitrogen, and pulverized using a cryoPREP tissue disruption system (CP02; Covaris). The pulverized tissue powder was divided into 10–20 mg aliquots for DNA, RNA, and protein extraction. Genomic DNA and RNA were extracted from cryopreserved tissue samples and blood samples using the AllPrep DNA/RNA Mini Kit, the QIAamp DNA Blood Mini Kit, the DNeasy Blood & Tissue Kit, and the RNeasy Micro Kit (Qiagen, Valencia, CA, USA) according to the manufacturers’ instructions. The concentration and integrity of the extracted RNA were assessed via a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Whole exome sequencing (WES) analysis To generate standard exome capture libraries, the SureSelect Human All Exon V6 + UTR probe set (Agilent Technologies) was used according to the manufacturer’s instructions. Briefly, genomic DNA extracted from the samples was fragmented and processed for adaptor ligation and PCR amplification. The exome capture libraries were constructed and amplified. The final purified products were quantified according to the qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms) and qualified using a TapeStation DNA Screentape D1000 (Agilent). The products were sequenced using the NovaSeq platform (Illumina Inc., San Diego, CA, USA). Quality-calibrated fastq files were aligned to the reference genome hg19 using BWA-MEM . Each read group was aligned to the reference genome separately, and all read group alignments that belonged to a single aliquot were merged using Picard Tools, SortSam and MergeSamFiles. Duplicate reads were flagged to prevent downstream variant call errors. Local realignment of insertions and deletions was performed using GATK:IndelRealigner . A base quality score recalibration (BQSR) step was then performed using GATK:BaseRecalibrator . Somatic mutation detection was performed via MuTect and GATK:Mutect2 , and the mutations were annotated with an oncotator . The GATK Haplotyper caller was used to detect germline variants , . Whole-transcriptome sequencing (WTS) Only high-quality RNA (RNA integrity ≥ 7.0) generated using the Illumina TruSeq Stranded mRNA Sample Prep Kit (Illumina) was used for RNA library construction. Qualified libraries were indexed, and paired-end 100-bp reads were sequenced using the Illumina NovaSeq platform. The quality of the sequencing reads was evaluated and examined for primer/adaptor sequence contamination via FastQC (v0.11.7) . Trim Galore was employed for trimming the reads, applying a threshold of an average sequence quality > 30. The trimmed reads were subsequently aligned to the human reference genome using STAR (STAR-2.7.0). Gene expression profiling was conducted with the HTseq-count algorithm, which uses the raw read counts as input data. Finally, edgeR was used to convert the count data to reads per kilobase million (RPKM), or the transcript per million (TPM) algorithm from the bioinfokit Python package (v1.0.4) was used to adjust different sequencing throughputs between the samples. Proteomic analysis Proteomic analysis was performed according to the standard protocol reported by the Clinical Proteomic Tumor Analysis Consortium (CPTAC) . Briefly, cryopreserved tissue powder samples were solubilized in sodium dodecyl sulfate (SDS) solubilization buffer (5% SDS, 50 mM TEAB; pH 8.5) using an S220 focused ultrasonicator (Covaris). Proteins were digested using S-Trap™ spin columns (ProtiFi, Huntington, NY, USA), cleaned with C18 spin columns (Thermo Fisher Scientific, Rockford, IL, USA), and the desalted peptide samples were labeled using TMT11plex reagents (Thermo Fisher Scientific). Pooled 11-plex tandem mass tag (TMT)-labeled samples were processed for peptide fractionation using an Agilent 1260 Infinity HPLC system (Agilent) and then prepared for global proteome and phosphoproteome analyses. Phosphopeptide enrichment was performed using immobilized metal affinity chromatography . The TMT-labeled peptides were loaded onto a trap column (Acclaum PepMap TM 100, 75 mm × 2 cm), separated on an analytical column (EASY-Spray column, 75 mm × 50 cm; Thermo Fisher Scientific) via an Ultimate 3000 RSLCnano system (Thermo Fisher Scientific), and analyzed using the top 10 data-dependent methods using a Q Exactive HF-X hybrid quadrupole–orbitrap mass spectrometer (Thermo Fisher Scientific). Cell lines and reagents Human breast cancer cell lines were cultured according to the distributors’ instructions. We obtained breast cancer cell lines from the American Type Culture Collection (ATCC) in Manassas, Virginia, USA; the Korean Cell Line Bank in Seoul, Korea; and the DSMZ in Braunschweig, Germany. MDA-MD-453, BT-474, and HCC-1954 cells were maintained in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS). SKBR3 cells were cultured in McCoy’s 5A medium supplemented with 10% FBS. JIMT-1 cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% FBS. Penicillin/streptomycin (1%) was added to all of the culture media, and the cells were cultured at 37 °C with 5% CO2. Droplet digital PCR (ddPCR) The copy number of the HER2 gene was analyzed using droplet digital PCR (ddPCR) on a QX200 Droplet Digital PCR System by Bio-Rad Laboratories. Fluorescent probes (FAM and HEX) were prepared from the PrimePCR TM ddPCR TM Copy Number Assay for ddPCR (dHsaCP1000116 for HER2 , dHsaCP2500349 for EIF2C1 , and dHsaCNS516206038 for POLR2A as the reference control) (Bio-Rad Laboratories Inc., Hercules, CA, USA). A total of 20 µl of PCR mixture was prepared with 10 ng of genomic DNA using 2X ddPCR supermix for the probe (10 µl) and 20X HER2 and EIF2C1 or POLR2A probes (FAM/HEX) (1 µl). PCRs were run on a Mastercycler nexus gradient Thermal Cycler (Eppendorf, Hamburg, Germany) at 95 °C for 10 min, followed by 40 cycles of 94 °C for 30 s, 55 °C for 60 s, and a 10 min incubation at 98 °C. Then, the PCR plates were read with a Bio-Rad QX200 droplet reader (Bio-Rad Laboratories, Inc.) and QuantaSoft TM version 1.4.0 software (Bio-Rad Laboratories, Inc.) to assess the number of droplets positive for HER2 and EIF2C1 or POLR2A . Viability assay Breast cancer cell lines were seeded in triplicate at 2–4 × 10 3 cells per well, depending on the cell line, in 96-well plates. We also treated breast cancer cell lines with neratinib (Selleckchem, Houston, TX, USA) or lapatinib (Selleckchem). After 48 h of treatment with the test compounds, cell viability was examined using a Cell Counting Kit-8 (CCK-8, Dojindo Laboratories, Kumamoto, Japan). Dimethyl sulfoxide (DMSO) (0.1%) was used as the control. The half-maximal inhibitory concentration (IC 50 ) was analyzed using GraphPad Prism (version 5.03, GraphPad Software, San Diego, CA, USA). Western blot Total protein was extracted from cell lines or tissues using RIPA lysis buffer (Thermo Fisher Scientific) with protease and phosphatase inhibitor cocktails (GenDEPOT, Baker, TX, USA) following the manufacturer’s instructions. The protein concentrations were determined via a protein assay solution from Bio-Rad Laboratories. Proteins were separated via sodium dodecyl sulfate‒polyacrylamide gel electrophoresis (SDS‐PAGE), and the separated proteins were transferred to polyvinylidene fluoride (PVDF) membranes (Merck Millipore, Billerica, MA, USA). The blots were then blocked for 1 h with 5% skim milk (BD Biosciences, San Diego, CA, USA) in Tris-buffered saline containing 0.1% Tween 20 and incubated overnight with the indicated primary antibodies. Antibodies against the following proteins were used: phospho-HER2-Tyr1248 (Cell Signaling Technology, CST), HER2 (Abcam, Cambridge, UK), phospho-HER2-Ser1054 (Thermo Fisher Scientific), and β-actin (CST, used as the loading control). Proteins were visualized via a horseradish peroxidase-conjugated secondary antibody (CST) and an enhanced chemiluminescence (ECL) reagent (Bio-Rad Laboratories). Multiplex immunofluorescence staining Multiplex immunofluorescence staining for cytokeratin (CK), PD-L1, CD68, CD8, FOXP-3, and PD-1 proteins was performed using a Leica Bond Rx Automated Stainer (Leica Biosystems, Newcastle, UK). Slides with tumor tissue sections were incubated with antibodies against the following proteins: CK (NBP2-29429; Novus Bio; dilution 1:300), PD-L1 (136845; CST; dilution 1:300), CD68 (ab192847; Abcam; dilution 1:300), CD8 (MCA1817; Bio-Rad; dilution 1:300), FOXP3 (ab20034; Abcam; dilution 1:100), and PD-1 (ab137132; Abcam; dilution 1:500). Then, the sections were incubated with the Polymer HRP Ms+Rb (ARH1001EA, AKOYA Biosciences) and subsequently treated with an Opal fluorophore (Opal690, Opal570, Opal520, Opal620, Opal480, or Opal780). Nuclei were subsequently visualized with DAPI, and the sections were coverslipped with ProLong Gold antifade reagent (P36934, Invitrogen). The slides were scanned using the Vectra Polaris Automated Quantitative Pathology Imaging System (Akoya Biosciences, Marlborough, MA, USA), and regions of interest were selected for all available tumor-containing areas. Images were analyzed using Inform 2.5 software (Akoya Biosciences, Marlborough, MA, USA) and TIBCO Spotfire™ (TIBCO, Palo Alto, CA, USA). Homologous recombination deficiency (HRD) score The HRD score was determined as a simple sum of three factors: NtAI (number of telomeric allelic imbalances), LST (large-scale state transitions), and HRD-LOH (HRD loss of heterozygosity). Tumors with HRD scores ≥ 42 were defined as HRD-high. For the determination of HRD scores via WES, BAM files of tumor samples were applied to the Sequenza followed by the scarHRD and/or HRDetect package using the default parameters as previously described – . Differentially Expressed Gene (DEG) analyses DEG analyses were conducted using edgeR , Cuffdiff2 , and DESeq2 following previously published protocols with default parameters. edgeR employed generalized linear models with tagwise dispersion, and raw counts were used as the primary input. As Cuffdiff2 does not work with count matrices, aligned transcriptomic reads were assembled into transcripts using Cufflinks with quartile normalization, bias correction, multiread correction, and a reference gene model. Cuffdiff2 identified DEGs via geometric library normalization and per-condition dispersion estimation. Adjusted p values were calculated with Benjamini‒Hochberg false discovery correction (5%) for all methods. Significance of mutation The MutSigCV algorithm (v.1.3.5) was employed to discern significantly mutated genes in our cohort. The default parameters of three inputs, coverage table files, covariate table files and mutation type dictionaries, were used to run MutSigCV, which uses a statistical framework to evaluate the observed mutation patterns against background mutation rates while considering covariates, including sequence context and replication timing. Integrative data analysis Subtypes were identified on the basis of individual types of data (mRNA sequencing data, whole proteome, phosphoproteome) expressed in at least 70% of patient samples. Median absolute deviations (MADs) were calculated for the individual data, and the top 10% of the most variably expressed genes were subjected to subsequent statistical analyses. On the basis of the MAD values, we then selected the molecules with the top 10% of MADs: 1,468 genes, 723 proteins, and 128 phosphopeptides. Next, we applied network integration based on consensus clustering to conduct unsupervised classification and to identify proteogenomic features that represent characteristic expression patterns for each cluster. The three data types were aggregated to generate an integrated network and similarity matrix, followed by the designation of consensus clusters. The criteria for determining the number of clusters were as follows: relatively high consistency within clusters, a relatively low variation coefficient, and no appreciable rise in the area under the cumulative distribution function (CDF) curve. The samples most representative of the clusters were identified on the basis of their positive silhouette width. To determine the representation of intrinsic subtypes of our cohort, the PAM50 model was applied to TPM-normalized RNA-seq data using the “genefu” R package (v3.12) . Single-sample gene set enrichment analysis (ssGSEA) was performed for each sample using log2-transformed and quantile-normalized RNA-seq data and the global proteome for the molecular signatures database (MSigDB) GeneSets via GenePattern 2.0 and Gene Ontology (GO) biological process signatures. The gene set scores represent the enrichment scores obtained from ssGSEA. To generate the receiver operating characteristic (ROC) curve and calculate the area under the curve (AUC) from the ssGSEA results and binary phenotype classification, we used ssGSEA_ROC from GenePattern 2.0. Statistical significance was determined using a false discovery rate (FDR) of less than 0.01. To obtain individual enrichment scores for each pairing of a sample and a gene set, we used the ssGSEA projection on the GenePattern server using t test statistics with a p- value less than 0.001. Immune profiling and downstream analysis To calculate RNA-based tumor immune scores and estimate immune-cell-specific contributions to each tumor, TPM data were analyzed using ESTIMATE , CIBERSORT in absolute mode , xCell , and the MCP counter . We also inferred the immune cell infiltration by ssGSEA using a published immune gene signature . Protein-based immune scores for stimulatory and inhibitory immune modulators, as well as the set of HLA proteins, were calculated as the mean of the protein log ratios in each set defined in Thorsson et al . Kinase phosphorylation outliers To identify the kinase activity characteristic of each PAM50 and major integrative clusters (iClusters), we used BlackSheep’s differential extremes value analysis module . For each phosphosite, the median and interquartile range (IQR) were calculated across all tumors. A site was defined as an outlier if it was more than 1.5 times the IQR above the median. The phosphosites were then collapsed into proteins by counting outlier and nonoutlier values per sample. For each group of interest (e.g., PAM50 or iCluster), proteins that were not enriched in outliers in that group and proteins without at least 30% of samples with an outlier were removed. Following filtering, outlier and nonoutlier sites per gene were counted for each group of interest, and Fisher’s exact test was used to calculate a p -value with multiple hypothesis correction via the Benjamini‒Hochberg (BH) procedure. For additional insight into kinase activity, we visualized the enrichment of kinase activation loop phosphorylation, which was calculated using a rank sum test. Kinase-substrate enrichment analysis (KSEA) KSEA was performed using the KSEA App web-based tool with phosphosite data with a cutoff of p < 0.05 and a substrate count of ≥ 1. Student’s t -test was used to calculate the p value, and the mean ratio was used to determine the fold change compared with the other four subsets. KSEA computes a normalized score to measure relative kinase activity in cluster 1 versus others using the difference in the mean log fold changes of predefined substrate groups of the kinase, as referenced from the PhosphoSitePlus and NetworKIN databases . In the latter database, kinase–substrate annotations with a NetworKIN score of less than five are discarded. The mean log fold changes of all substrates were divided by the standard deviation of log fold changes across all phosphosites, followed by an adjustment to the score for the number of substrates in the data that were annotated to the kinase. The statistical significance of the scores was calculated using a one-tailed z test, followed by BH multiple testing correction. Kinases with z scores and corresponding p values less than 0.05 were considered significant. 

A total of 178 participants with histologically defined breast cancer, aged 40 years or younger, treated at the National Cancer Center (NCC) in Korea were included in this study. Tumor and adjacent normal tissue samples were obtained from surgically resected specimens, and blood samples were collected from the patients. This study was approved by the Institutional Review Board of the NCC (IRB nos. NCCNCS 13717, NCC2017-0127, and NCC2020-0135). Since some patients had limited samples, only 126 patients for whom all proteogenomic data were fully generated were used for analysis. We retrospectively reviewed the medical and pathological records of the patients, including histological diagnoses of surgical specimens, tumor staging, histological grade, treatment history (type of surgery, use of chemotherapy, hormone therapy, anti-HER2 therapy, and radiotherapy), and follow-up data (recurrence, metastasis, and death). Tissue specimens and blood samples were provided by the NCC Bio Bank of the National Cancer Center, Korea. The frozen tissue samples were weighed and washed with cold phosphate-buffered saline (PBS) to remove blood contamination, placed in a tube (Covaris, Woburn, MA, USA), snap-frozen in liquid nitrogen, and pulverized using a cryoPREP tissue disruption system (CP02; Covaris). The pulverized tissue powder was divided into 10–20 mg aliquots for DNA, RNA, and protein extraction. Genomic DNA and RNA were extracted from cryopreserved tissue samples and blood samples using the AllPrep DNA/RNA Mini Kit, the QIAamp DNA Blood Mini Kit, the DNeasy Blood & Tissue Kit, and the RNeasy Micro Kit (Qiagen, Valencia, CA, USA) according to the manufacturers’ instructions. The concentration and integrity of the extracted RNA were assessed via a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). 

To generate standard exome capture libraries, the SureSelect Human All Exon V6 + UTR probe set (Agilent Technologies) was used according to the manufacturer’s instructions. Briefly, genomic DNA extracted from the samples was fragmented and processed for adaptor ligation and PCR amplification. The exome capture libraries were constructed and amplified. The final purified products were quantified according to the qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms) and qualified using a TapeStation DNA Screentape D1000 (Agilent). The products were sequenced using the NovaSeq platform (Illumina Inc., San Diego, CA, USA). Quality-calibrated fastq files were aligned to the reference genome hg19 using BWA-MEM . Each read group was aligned to the reference genome separately, and all read group alignments that belonged to a single aliquot were merged using Picard Tools, SortSam and MergeSamFiles. Duplicate reads were flagged to prevent downstream variant call errors. Local realignment of insertions and deletions was performed using GATK:IndelRealigner . A base quality score recalibration (BQSR) step was then performed using GATK:BaseRecalibrator . Somatic mutation detection was performed via MuTect and GATK:Mutect2 , and the mutations were annotated with an oncotator . The GATK Haplotyper caller was used to detect germline variants , . 

Only high-quality RNA (RNA integrity ≥ 7.0) generated using the Illumina TruSeq Stranded mRNA Sample Prep Kit (Illumina) was used for RNA library construction. Qualified libraries were indexed, and paired-end 100-bp reads were sequenced using the Illumina NovaSeq platform. The quality of the sequencing reads was evaluated and examined for primer/adaptor sequence contamination via FastQC (v0.11.7) . Trim Galore was employed for trimming the reads, applying a threshold of an average sequence quality > 30. The trimmed reads were subsequently aligned to the human reference genome using STAR (STAR-2.7.0). Gene expression profiling was conducted with the HTseq-count algorithm, which uses the raw read counts as input data. Finally, edgeR was used to convert the count data to reads per kilobase million (RPKM), or the transcript per million (TPM) algorithm from the bioinfokit Python package (v1.0.4) was used to adjust different sequencing throughputs between the samples. 

Proteomic analysis was performed according to the standard protocol reported by the Clinical Proteomic Tumor Analysis Consortium (CPTAC) . Briefly, cryopreserved tissue powder samples were solubilized in sodium dodecyl sulfate (SDS) solubilization buffer (5% SDS, 50 mM TEAB; pH 8.5) using an S220 focused ultrasonicator (Covaris). Proteins were digested using S-Trap™ spin columns (ProtiFi, Huntington, NY, USA), cleaned with C18 spin columns (Thermo Fisher Scientific, Rockford, IL, USA), and the desalted peptide samples were labeled using TMT11plex reagents (Thermo Fisher Scientific). Pooled 11-plex tandem mass tag (TMT)-labeled samples were processed for peptide fractionation using an Agilent 1260 Infinity HPLC system (Agilent) and then prepared for global proteome and phosphoproteome analyses. Phosphopeptide enrichment was performed using immobilized metal affinity chromatography . The TMT-labeled peptides were loaded onto a trap column (Acclaum PepMap TM 100, 75 mm × 2 cm), separated on an analytical column (EASY-Spray column, 75 mm × 50 cm; Thermo Fisher Scientific) via an Ultimate 3000 RSLCnano system (Thermo Fisher Scientific), and analyzed using the top 10 data-dependent methods using a Q Exactive HF-X hybrid quadrupole–orbitrap mass spectrometer (Thermo Fisher Scientific). 

Human breast cancer cell lines were cultured according to the distributors’ instructions. We obtained breast cancer cell lines from the American Type Culture Collection (ATCC) in Manassas, Virginia, USA; the Korean Cell Line Bank in Seoul, Korea; and the DSMZ in Braunschweig, Germany. MDA-MD-453, BT-474, and HCC-1954 cells were maintained in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS). SKBR3 cells were cultured in McCoy’s 5A medium supplemented with 10% FBS. JIMT-1 cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% FBS. Penicillin/streptomycin (1%) was added to all of the culture media, and the cells were cultured at 37 °C with 5% CO2. 

The copy number of the HER2 gene was analyzed using droplet digital PCR (ddPCR) on a QX200 Droplet Digital PCR System by Bio-Rad Laboratories. Fluorescent probes (FAM and HEX) were prepared from the PrimePCR TM ddPCR TM Copy Number Assay for ddPCR (dHsaCP1000116 for HER2 , dHsaCP2500349 for EIF2C1 , and dHsaCNS516206038 for POLR2A as the reference control) (Bio-Rad Laboratories Inc., Hercules, CA, USA). A total of 20 µl of PCR mixture was prepared with 10 ng of genomic DNA using 2X ddPCR supermix for the probe (10 µl) and 20X HER2 and EIF2C1 or POLR2A probes (FAM/HEX) (1 µl). PCRs were run on a Mastercycler nexus gradient Thermal Cycler (Eppendorf, Hamburg, Germany) at 95 °C for 10 min, followed by 40 cycles of 94 °C for 30 s, 55 °C for 60 s, and a 10 min incubation at 98 °C. Then, the PCR plates were read with a Bio-Rad QX200 droplet reader (Bio-Rad Laboratories, Inc.) and QuantaSoft TM version 1.4.0 software (Bio-Rad Laboratories, Inc.) to assess the number of droplets positive for HER2 and EIF2C1 or POLR2A . 

Breast cancer cell lines were seeded in triplicate at 2–4 × 10 3 cells per well, depending on the cell line, in 96-well plates. We also treated breast cancer cell lines with neratinib (Selleckchem, Houston, TX, USA) or lapatinib (Selleckchem). After 48 h of treatment with the test compounds, cell viability was examined using a Cell Counting Kit-8 (CCK-8, Dojindo Laboratories, Kumamoto, Japan). Dimethyl sulfoxide (DMSO) (0.1%) was used as the control. The half-maximal inhibitory concentration (IC 50 ) was analyzed using GraphPad Prism (version 5.03, GraphPad Software, San Diego, CA, USA). 

Total protein was extracted from cell lines or tissues using RIPA lysis buffer (Thermo Fisher Scientific) with protease and phosphatase inhibitor cocktails (GenDEPOT, Baker, TX, USA) following the manufacturer’s instructions. The protein concentrations were determined via a protein assay solution from Bio-Rad Laboratories. Proteins were separated via sodium dodecyl sulfate‒polyacrylamide gel electrophoresis (SDS‐PAGE), and the separated proteins were transferred to polyvinylidene fluoride (PVDF) membranes (Merck Millipore, Billerica, MA, USA). The blots were then blocked for 1 h with 5% skim milk (BD Biosciences, San Diego, CA, USA) in Tris-buffered saline containing 0.1% Tween 20 and incubated overnight with the indicated primary antibodies. Antibodies against the following proteins were used: phospho-HER2-Tyr1248 (Cell Signaling Technology, CST), HER2 (Abcam, Cambridge, UK), phospho-HER2-Ser1054 (Thermo Fisher Scientific), and β-actin (CST, used as the loading control). Proteins were visualized via a horseradish peroxidase-conjugated secondary antibody (CST) and an enhanced chemiluminescence (ECL) reagent (Bio-Rad Laboratories). 

Multiplex immunofluorescence staining for cytokeratin (CK), PD-L1, CD68, CD8, FOXP-3, and PD-1 proteins was performed using a Leica Bond Rx Automated Stainer (Leica Biosystems, Newcastle, UK). Slides with tumor tissue sections were incubated with antibodies against the following proteins: CK (NBP2-29429; Novus Bio; dilution 1:300), PD-L1 (136845; CST; dilution 1:300), CD68 (ab192847; Abcam; dilution 1:300), CD8 (MCA1817; Bio-Rad; dilution 1:300), FOXP3 (ab20034; Abcam; dilution 1:100), and PD-1 (ab137132; Abcam; dilution 1:500). Then, the sections were incubated with the Polymer HRP Ms+Rb (ARH1001EA, AKOYA Biosciences) and subsequently treated with an Opal fluorophore (Opal690, Opal570, Opal520, Opal620, Opal480, or Opal780). Nuclei were subsequently visualized with DAPI, and the sections were coverslipped with ProLong Gold antifade reagent (P36934, Invitrogen). The slides were scanned using the Vectra Polaris Automated Quantitative Pathology Imaging System (Akoya Biosciences, Marlborough, MA, USA), and regions of interest were selected for all available tumor-containing areas. Images were analyzed using Inform 2.5 software (Akoya Biosciences, Marlborough, MA, USA) and TIBCO Spotfire™ (TIBCO, Palo Alto, CA, USA). 

The HRD score was determined as a simple sum of three factors: NtAI (number of telomeric allelic imbalances), LST (large-scale state transitions), and HRD-LOH (HRD loss of heterozygosity). Tumors with HRD scores ≥ 42 were defined as HRD-high. For the determination of HRD scores via WES, BAM files of tumor samples were applied to the Sequenza followed by the scarHRD and/or HRDetect package using the default parameters as previously described – . 

DEG analyses were conducted using edgeR , Cuffdiff2 , and DESeq2 following previously published protocols with default parameters. edgeR employed generalized linear models with tagwise dispersion, and raw counts were used as the primary input. As Cuffdiff2 does not work with count matrices, aligned transcriptomic reads were assembled into transcripts using Cufflinks with quartile normalization, bias correction, multiread correction, and a reference gene model. Cuffdiff2 identified DEGs via geometric library normalization and per-condition dispersion estimation. Adjusted p values were calculated with Benjamini‒Hochberg false discovery correction (5%) for all methods. 

The MutSigCV algorithm (v.1.3.5) was employed to discern significantly mutated genes in our cohort. The default parameters of three inputs, coverage table files, covariate table files and mutation type dictionaries, were used to run MutSigCV, which uses a statistical framework to evaluate the observed mutation patterns against background mutation rates while considering covariates, including sequence context and replication timing. 

Subtypes were identified on the basis of individual types of data (mRNA sequencing data, whole proteome, phosphoproteome) expressed in at least 70% of patient samples. Median absolute deviations (MADs) were calculated for the individual data, and the top 10% of the most variably expressed genes were subjected to subsequent statistical analyses. On the basis of the MAD values, we then selected the molecules with the top 10% of MADs: 1,468 genes, 723 proteins, and 128 phosphopeptides. Next, we applied network integration based on consensus clustering to conduct unsupervised classification and to identify proteogenomic features that represent characteristic expression patterns for each cluster. The three data types were aggregated to generate an integrated network and similarity matrix, followed by the designation of consensus clusters. The criteria for determining the number of clusters were as follows: relatively high consistency within clusters, a relatively low variation coefficient, and no appreciable rise in the area under the cumulative distribution function (CDF) curve. The samples most representative of the clusters were identified on the basis of their positive silhouette width. To determine the representation of intrinsic subtypes of our cohort, the PAM50 model was applied to TPM-normalized RNA-seq data using the “genefu” R package (v3.12) . Single-sample gene set enrichment analysis (ssGSEA) was performed for each sample using log2-transformed and quantile-normalized RNA-seq data and the global proteome for the molecular signatures database (MSigDB) GeneSets via GenePattern 2.0 and Gene Ontology (GO) biological process signatures. The gene set scores represent the enrichment scores obtained from ssGSEA. To generate the receiver operating characteristic (ROC) curve and calculate the area under the curve (AUC) from the ssGSEA results and binary phenotype classification, we used ssGSEA_ROC from GenePattern 2.0. Statistical significance was determined using a false discovery rate (FDR) of less than 0.01. To obtain individual enrichment scores for each pairing of a sample and a gene set, we used the ssGSEA projection on the GenePattern server using t test statistics with a p- value less than 0.001. 

To calculate RNA-based tumor immune scores and estimate immune-cell-specific contributions to each tumor, TPM data were analyzed using ESTIMATE , CIBERSORT in absolute mode , xCell , and the MCP counter . We also inferred the immune cell infiltration by ssGSEA using a published immune gene signature . Protein-based immune scores for stimulatory and inhibitory immune modulators, as well as the set of HLA proteins, were calculated as the mean of the protein log ratios in each set defined in Thorsson et al . 

To identify the kinase activity characteristic of each PAM50 and major integrative clusters (iClusters), we used BlackSheep’s differential extremes value analysis module . For each phosphosite, the median and interquartile range (IQR) were calculated across all tumors. A site was defined as an outlier if it was more than 1.5 times the IQR above the median. The phosphosites were then collapsed into proteins by counting outlier and nonoutlier values per sample. For each group of interest (e.g., PAM50 or iCluster), proteins that were not enriched in outliers in that group and proteins without at least 30% of samples with an outlier were removed. Following filtering, outlier and nonoutlier sites per gene were counted for each group of interest, and Fisher’s exact test was used to calculate a p -value with multiple hypothesis correction via the Benjamini‒Hochberg (BH) procedure. For additional insight into kinase activity, we visualized the enrichment of kinase activation loop phosphorylation, which was calculated using a rank sum test. 

KSEA was performed using the KSEA App web-based tool with phosphosite data with a cutoff of p < 0.05 and a substrate count of ≥ 1. Student’s t -test was used to calculate the p value, and the mean ratio was used to determine the fold change compared with the other four subsets. KSEA computes a normalized score to measure relative kinase activity in cluster 1 versus others using the difference in the mean log fold changes of predefined substrate groups of the kinase, as referenced from the PhosphoSitePlus and NetworKIN databases . In the latter database, kinase–substrate annotations with a NetworKIN score of less than five are discarded. The mean log fold changes of all substrates were divided by the standard deviation of log fold changes across all phosphosites, followed by an adjustment to the score for the number of substrates in the data that were annotated to the kinase. The statistical significance of the scores was calculated using a one-tailed z test, followed by BH multiple testing correction. Kinases with z scores and corresponding p values less than 0.05 were considered significant. 

Patient Cohort Description The baseline characteristics of the 126 patients with YBC are summarized in Supplementary Table . The median patient age was 37 years (range: 25–40). Seventy-six patients (60.3%) had stage II disease, and 64 patients (50.8%) had node metastasis at the time of diagnosis. Sixty-five patients (51.6%) had a histological grade of 3. ER was positive in 84 patients (66.7%), PR was positive in 73 patients (57.9%), and HER2 was positive in 25 patients (19.8%). Based on the clinical classification of breast cancer with ER, PR, HER2, and Ki-67 statuses, 33 patients were classified as luminal A type (26.2%), 55 as luminal B type (43.7%), 10 as HER2 type (7.9%), and 25 as TNBC type (19.8%). All patients were diagnosed with invasive ductal carcinoma. Overall, 113 patients (89.7%) underwent systemic chemotherapy, 104 patients (82.5%) underwent adjuvant radiotherapy, and 89 patients (70.6%) underwent adjuvant hormone therapy. The median follow-up period was 94.3 months (range: 4.0–254.3), 37 patients (29.4%) experienced recurrence, and 17 (13.5%) died. All of the patients in this study ( n = 126, 100%) were initially diagnosed under the age of 40 and premenopausal, and TCGA data largely consists of patients with breast cancer over the age of 40 ( n = 1,998, 91.1%). Compared with the TCGA cohort ( n = 1116), our patients demonstrated a greater proportion of the basal subtype (28.6% KNCC-YBC vs. 16.5% TCGA) and HER2-enriched subtype (15.9% vs. 6.3%). This difference was maintained when compared to young patients (age ≤ 40) in the TCGA cohort (basal type 17.7%, HER2 type 4.4% in TCGA_YBC [ n = 158]). This finding is comparable to the findings of a previous study in which Asian patients with breast cancer had a greater prevalence of HER2-positive cancers than Caucasian women did , . Genomic landscape of early-onset breast cancer We performed an integrated analysis of molecular data from patients through WES, WTS, and mass spectrometry-based proteomics, which provided robust quantification of proteins and phosphoproteins (Fig. ). Whole-exome sequencing of 126 samples revealed 30,193 somatic mutations in the study cohort. Each patient harbored a median of 119 somatic mutations per sample (range: 2–2858 somatic mutations). Significance analysis of the somatic mutations, employing covariates related to the mutational processes, identified the key mutated genes, including TP53, PIK3CA, GATA3, MAP2K4 , and ARID1A (Fig. ). TP53 (60%), PIK3CA (43%), and GATA3 (36%) were the most frequently mutated genes in this cohort. GATA3 somatic mutations were mostly mutually exclusive with TP53 mutations (Fisher’s exact test, P < 0.01), with PIK3CA mutations exhibiting similar mutual exclusivity patterns ( P < 0.01). Overall, approximately 68% of the patients possessed at least one mutation in one of these three genes. The prevalence of somatic mutations in MAP2K4 and ARID1A in breast cancer has also been documented in other studies – , indicating the presence of a congruent mutational profile within the early-onset Asian breast cancer demographic. Molecular subtype classification on the basis of multiomic profiles In addition to classical pathological and clinical laboratory-based stratification of the molecular subtypes of the cohorts using hormone receptor and HER2 immunohistochemistry (IHC) analyses, we attempted a robust molecular classification of YBC cohorts based on a multiomics dataset. Unsupervised clustering of the three data platforms (proteome, transcriptome, and phosphoproteome) resulted in discordant stratification unique to each molecular layer. A consensus integrative classification encompassing the three molecular data layers produced five iClusters within the young Korean breast cancer cohort (Fig. ). Cluster 1 was significantly concordant with the classical PAM50-based HER2-enriched subtype. This cluster was characterized by a high prevalence of HER2 amplification and TP53 mutations (69.6% and 82.6% in cluster 1, respectively). The key pathways in this cluster included protein secretion, DNA replication, and repair. Cluster 2 was almost invariably consistent with the basal-like cancer type and presented a high frequency of TP53 mutations and an absence of hormone receptor and HER2 mutations. Cluster 2 was associated with the activation of pathways such as KRAS signaling, autophagy, and cell cycle processes and the negative regulation of apoptosis. Clusters 1 and 2 demonstrated the worst prognosis, with recurrence and metastasis events frequently observed within 5 years after surgical resection (Fig. ). Cluster 3 largely overlapped with the classical luminal B subtype, whereas cluster 4 demonstrated closer alignment with the classical luminal A subtype. The important pathways of cluster 3 and cluster 4 included estrogen signaling, the cell cycle, DNA replication and repair, and fatty acid metabolism (Fig. ). Those clusters of young Korean patients with breast cancer experienced more frequent recurrence and/or death events 5 or more years post-diagnosis of breast cancer (Fig. ). All cases of late relapse in these clusters were ipsilateral (11/11, P < 0.00001, Fisher’s exact test); thus, the late recurrent tumors seemed to have originated from residual and/or resistant tumor cells left from their primary tumors after a sustained period of dormancy rather than from the manifestation of a second primary tumor. In contrast, basal-type and HER2-type Korean YBC demonstrated completely opposite biphasic patterns of survival probability; these groups’ relapse and consequent death events were more clustered in the first 5 years, followed by a plateau in survival probability. The dichotomy of biphasic survival probability between luminal and nonluminal Korean patients with breast cancer is substantially distinct from the clinical course of Caucasian breast cancer, where luminal A breast cancer shows the most indolent phenotype, followed by luminal B breast cancer and HER2 breast cancer, with the basal group having the highest risk patterns . The results of the integrated analysis revealed a distinct subtype, cluster 5, which had a very good prognosis, unlike the other four clusters (Fig. ). The significant pathways of cluster 5 included vasculature development, endothelial cell proliferation, and oxidative phosphorylation. To examine the impact of somatic copy number alterations (SCNAs) on the levels of mRNAs and proteins, we conducted univariate correlation analysis to identify features that exhibited statistically significant correlations (Fig. ). This analysis considered changes in ‘cis’—affecting genes within the altered loci—and ‘trans’—impacting genes located elsewhere. We evaluated the genes in this cohort, each with corresponding SCNA, mRNA, and protein data, by calculating the correlation and its statistical significance for all CNA–mRNA and CNA–protein pairings. Notably, 61% of the CNA-mRNA pairs and 30% of the CNA–protein pairs demonstrated significant positive ‘cis’ correlations. We further scrutinized the subset of these pairs that included well-characterized oncogenes and tumor suppressor genes. ‘Trans’ effects appeared as prominent vertical bands in the analysis, with frequency histograms highlighting significant ‘hot spots’ of ‘trans’ correlation. Interestingly, correlations between CNA and proteins were generally weaker than those observed for CNA and mRNAs. The most pronounced ‘trans’ relationships at the protein level were identified in specific CNA regions, including regions of chromosome 5 and chromosome 8, in early-onset Korean breast cancer. Targetable elements associated with clusters and genomic alterations Next, we examined the functional impact of genetic alterations on the abundance of both cis -acting (cognate gene products) and trans- acting (other gene products) proteins and their phosphorylation sites (Fig. ). Several cis and trans associations were discernible. These include enhanced expression of HER2 proteins (Fig. ) and strong cis -phosphorylation (Fig. ) in HER2-positive breast cancers, cementing their role as potent drivers of alterations in a subset of breast cancers. Migration and invasion enhancer I ( MIEN1 ) and growth factor receptor-bound protein 7 ( GRB7 ) are genes located proximal to the ERBB2 gene and are concomitantly overexpressed in HER2-positive breast cancer subgroups (Fig. ). MIEN1 has been shown to functionally increase the invasive and migratory phenotypes of various types of cells in breast cancer . Additionally, GRB7 overexpression facilitates the phosphorylation of both AKT and HER2/neu in HER2/neu-overexpressing cells, promoting tumor growth in vivo . These findings suggest additional roles for these proximal genes in breast cancer. Several phosphorylation sites associated with common breast cancer oncogenes were identified and clustered in the cell cycle and DNA repair activities, which is consistent with their established roles in cancer. We investigated which kinase activities were enriched in the molecular subgroups defined in this study (Fig. ). As expected, HER2 kinase activity was upregulated in cluster 1, which was highly correlated with the classical PAM50 Her2 subgroup. Cluster 2 (basal subgroup) was enriched with cell cycle kinases, including CDC7, CDK1, and CDK2. Hormone receptor-positive cancer groups (clusters 3-4) presented increased activity of kinases involved in immune regulation. To identify potential therapeutic targets specific to each iCluster subtype, we used phosphoproteomic data to determine kinase activation (Fig. ). We identified phosphorylated kinases that were enriched in each iCluster subtype via outlier enrichment analysis. The dataset revealed several druggable targets and kinases, including ERBB2 (trastuzumab), in the cluster 1 subtype. The cluster 2 subtype includes isolated PML (tretinoin) and RCSD1 (imatinib). PDCD4 (paclitaxel) was found in cluster 4, whereas CD34 (prednisolone) and TJP1 (dexamethasone) were present in cluster 5. These proteins and kinases are potential targets for subtype-specific treatment. We also investigated the associations between specific genomic alterations and kinase activity (Fig. ). We identified 16 kinases significantly associated with genomic alterations. RB1, PTEN , and MAP2K4 mutations formed a strong subcluster with increased activity of several canonical cell cycle kinases, including the CDK1, CDK2, CDK7, and CDC7 proteins. These associations were also detected in ARID1A- mutant breast cancers. The PLK4, CSNK2A2, HER2, PRKCT, and PRKCH proteins are associated with HER2 genomic alterations. GATA3 genomic alteration was associated with increased activity of the HCK and TYK2 kinases, whereas PIK3CA genomic alteration was linked with the CDK19 kinase. Proteogenomic analysis of HER2-type YBC We assessed the HER2 driver status by integrating genomic and proteomic data to refine its classification and clinical significance. HER2 status was also grouped by classical IHC and in situ hybridization testing following ASCO-CAP guidelines (HER2+: IHC score of 3+ or an IHC score of 2+ fluorescence in situ hybridization (FISH); silver-enhanced in situ hybridization (SISH): HER2 equivocal status: IHC score of 2+ without FISH results or amplification by FISH without IHC results) . Proteogenomic analysis revealed that 19 patients who presented high levels of both HER2 protein and phosphorylated HER2 were designated as HER2 proteogenomic status positive (HER2 PG+) (Fig. ). HER2 PG+ cases were enriched mostly in cluster 1, followed by cluster 3. The majority of HER2 PG+ samples (16 out of 19) demonstrated largely consistent pathological designations with confirmed HER2-positive expression according to IHC, FISH, and SISH tests. The remaining three patients who were grouped according to HER2 equivocal status demonstrated a high HER2 copy number as a result of droplet digital PCR analysis. Although we identified 28 patients with HER2 gene amplification on the basis of copy number alterations, only 19 of these patients were positive for HER2 by PGstatus. Consistent with the proximal location of the HER2 gene on chromosome 17q12, the copy number, RNA, and protein expression of GRB7 and STARD3 were concomitantly increased in patients with HER2 amplification (Fig. ). The protein levels and phosphorylation metrics of HER2 more robustly and clearly segregated the designation between patients with positive (HER2-PG+) and negative (HER2-PG-) HER2 PG status (Fig. ). In the stratification of patients with breast cancer with HER2 driver dependency, there was a notable incompatibility between conventional PAM50-based classification, classical pathological HER2 type designation, and proteogenomic-based integrative clustering. To gain deeper insights into the biological characteristics that led to sample clustering within the HER2-enriched subtype (HER2E) group despite inconsistent ERBB2 status, we conducted an analysis of phosphosites that are known to canonically relay kinase signaling from the human ERBB2 protein. Patients assigned to the HER2E group were divided into HER2EPG- and HER2EPG+ groups on the basis of their PG status (Fig. ). As anticipated, all PAM50 HER2EPG+ samples presented elevated levels of ERBB2 phosphopeptides. Five additional cases of the PAM50 HER2-enriched subtype were classified into the HER2EPG- subgroup. A close examination of the five cases revealed that these samples shared common activated receptor tyrosine kinase (RTK) downstream pathway signatures, including ARAF, MAPK3, and AKT1 (Fig. ). This finding indicates that the PAM50 simplification scheme mistakenly designates some samples with activated RTKs as belonging to the HER2 subtype, regardless of their subtype-specific nature. Indeed, some samples demonstrated an activated phosphorylation profile of IGF2R, further suggesting that some breast cancer samples might be targetable with precise segregation of RTKs with a suitable driver dependency nature . Quantitative phosphoproteomic analysis revealed that HER2 PG+ cases exhibited elevated levels of HER2 serine phosphopeptides, including serine residues 1054 and 1066, in comparison to PG- cases (Fig. ). Given the close association between the phosphorylation status at serine 1054 (S-1054) and HER2 PG status, we investigated its potential as a predictive biomarker for the sensitivity of breast cancer cell lines to HER2 inhibitors. Despite the limited number of cell lines analyzed, the phosphorylation of S-1054 was found to be correlated with the response to lapatinib and neratinib treatments (Fig. ). We then compared the effects of a HER2 inhibitor on the phosphorylation status of HER2-positive breast cancer cell lines. The phosphorylation level of HER2 at tyrosine 1248 (Y-1248) decreased following treatment with neratinib. The phosphorylation of S-1054 decreased in HER2-positive breast cancer cell lines following neratinib treatment (Supplementary Fig. ). These results suggest that decreased phosphorylation of S-1054 on HER2 could serve as a predictive marker for sensitivity to HER2 inhibitors. Among patients with amplified HER2, the proteogenomic status of HER2 was potentially prognostic, whereas PAM50 HER2 enrichment was not (Fig. ). The patients with HER2 PG+ status tended to have poorer clinical outcomes among the patients with amplified HER2 ( p = 0.024). Proteogenomic analysis of the tumor immune microenvironment Using transcriptome, proteome, and phosphoproteome datasets, we attempted to classify patients according to their immune signals to understand their tumor immune microenvironments at the molecular and cellular levels (Fig. ). Clusters 1 and 2 were enriched in immune signatures, and the mRNA and protein expression levels of immune checkpoint targets, such as CD276, TIGIT, LAG3, and CTLA4, were also elevated relative to those in the other clusters. Cluster 3 included an increased proportion of patients with fewer immunogenic features, with decreased numbers of T, B, NK, and monocytic lineages, as well as immune checkpoint targets. Multiplex IHC staining of five different immune cell markers (CD8, CD68, FOXP3, PD1, and PD-L1) confirmed the validity of inferring tumor immune microenvironment features from a multiomics dataset (Fig. ). Correlation analyses of multiplexed IHC staining data revealed that YBC immune cells have highly intercorrelated molecular features, confirming the established correlative nature of immune checkpoint molecules, such as augmented PD-L1 staining, which is highly correlated with an increased number of exhausted CD8+ Tregs and immunosuppressive macrophages in the tumor parenchyma (Supplementary Fig. ). CD8-positive T-cell IHC further validated the inferred active immunogenic tumor microenvironment (Fig. ). We also identified three immune-based subtypes with distinct immune characteristics (Fig. ). Among the three clusters, immune-cluster 2 was assigned an immune hot phenotype, immune-cluster 3 an immune cold phenotype, and immune-cluster 1 an immune -intermediate phenotype. There was a positive correlation between the stromal and immune scores according to the three immune-based subtypes (R 2 = 0.4246, p < 0.0001; Fig. ). Classical RNA-seq-based PAM50 classification demonstrated that each subtype of YBC comprises a subset of patients with an active immune signature (Supplementary Fig. ), albeit with higher immune scores in the basal PAM50 subtype (Supplementary Fig. ) (ANOVA, P = 0.007). The active immune tumor microenvironment consisted of elevated levels of stimulatory immune proteins and HLA proteins and increased proportions of tumor-infiltrating lymphocytes (TILs) and myeloid dendritic cells (Supplementary Fig. ). In contrast to the PAM50 stratification, integrative multiomics-based clustering generated immune microenvironment features that segregated cluster 3 from the other clusters (ANOVA, P = 0.001) (Supplementary Fig. ). The evaluation of the clinical characteristics of the poorly immunogenic subgroup tended to be associated with a poor clinical prognosis (Supplementary Fig. ). The clinical association of poor immunogenic features with a poor prognosis was observed in the luminal B PAM50 subtypes in this study cohort (Fig. ), suggesting that immune microenvironment features in specific subgroups of early-onset breast cancer identified through multiomics integration may serve as prognostic molecular markers for stratifying high-risk late-onset luminal B subtype patients. Proteogenomic analysis of homologous recombination-deficient (HRD) breast cancer To investigate the mutational processes involved in early-onset East Asian breast cancer, we deconvoluted the mutational signatures in the cohort. Four prevalent mutational signatures—BRCAness, APOBEC, mismatch repair defects, and aging signatures—were observed (Fig. ). We also examined germline pathogenic variations, defined as mutations that were reported as the most likely disease-causing variants in ClinVar or mutations that are truncating proteins in 10 genes reported to increase breast cancer susceptibility with a relatively strong penetrance , , . Among the 10 breast cancer susceptibility genes, 10.3% (13/126) of patients harbored pathogenic germline variations, whereas one patient was observed to possess pathogenic variations in both RAD50 and BRCA1. Pathogenic BRCA1/2 germline mutations were more prevalent in Korean patients with YBC, affecting 8.7% of patients in the study cohort compared with 3.8% of patients in the TCGA database (Fisher’s test, P < 0.02). BRCA germline mutations, which are classical biomarkers of BRCAness and PARP inhibitors, only partially explained the proportion of patients with elevated genomic HRD levels (Fig. ). We interrogated the proteome data to search for proteins closely associated with the features of HRD (Fig. ). Among these, 20 highly ranked proteins were chosen whose combination could segregate genomic HRD-positive and HRD-negative cancers (Fig. ). Prognostic significance of protein‒RNA correlations Positive correlations existed between the majority of the quantified proteins (>90%) and their corresponding mRNA transcripts within the present study cohort, with a median of 0.34 (Supplementary Fig. ). These correlations were not contingent upon technical factors such as the average protein precursor area or the number of peptide spectral matches. Structural ribosomal proteins and those involved in translational processes presented the lowest correlations between protein and transcript abundance (Supplementary Fig. ), highlighting their function and regulatory control, primarily at the protein level. Notably, proteins derived from transcripts identified in breast cancer drivers, including estrogen receptor responses, exhibited the strongest mRNA‒protein correlations (Supplementary Fig. ), underscoring a direct connection to the protein phenotype in the major canonical breast cancer driver signaling pathway. In hormone receptor-positive breast cancer patients, the subgroup with a high RNA‒protein correlation was associated with a poorer prognosis than the subgroup with a low RNA‒protein correlation (Supplementary Fig. , p = 0.0049). This is reinforced by the greater correlation of genes that respond to estrogen receptor signaling blockade (Supplementary Fig. ). These findings suggest that protein‒RNA correlations could serve as promising prognostic markers for predicting late recurrence in the luminal subtype, necessitating comprehensive external validation accounting for potential confounders in future investigations. 

The baseline characteristics of the 126 patients with YBC are summarized in Supplementary Table . The median patient age was 37 years (range: 25–40). Seventy-six patients (60.3%) had stage II disease, and 64 patients (50.8%) had node metastasis at the time of diagnosis. Sixty-five patients (51.6%) had a histological grade of 3. ER was positive in 84 patients (66.7%), PR was positive in 73 patients (57.9%), and HER2 was positive in 25 patients (19.8%). Based on the clinical classification of breast cancer with ER, PR, HER2, and Ki-67 statuses, 33 patients were classified as luminal A type (26.2%), 55 as luminal B type (43.7%), 10 as HER2 type (7.9%), and 25 as TNBC type (19.8%). All patients were diagnosed with invasive ductal carcinoma. Overall, 113 patients (89.7%) underwent systemic chemotherapy, 104 patients (82.5%) underwent adjuvant radiotherapy, and 89 patients (70.6%) underwent adjuvant hormone therapy. The median follow-up period was 94.3 months (range: 4.0–254.3), 37 patients (29.4%) experienced recurrence, and 17 (13.5%) died. All of the patients in this study ( n = 126, 100%) were initially diagnosed under the age of 40 and premenopausal, and TCGA data largely consists of patients with breast cancer over the age of 40 ( n = 1,998, 91.1%). Compared with the TCGA cohort ( n = 1116), our patients demonstrated a greater proportion of the basal subtype (28.6% KNCC-YBC vs. 16.5% TCGA) and HER2-enriched subtype (15.9% vs. 6.3%). This difference was maintained when compared to young patients (age ≤ 40) in the TCGA cohort (basal type 17.7%, HER2 type 4.4% in TCGA_YBC [ n = 158]). This finding is comparable to the findings of a previous study in which Asian patients with breast cancer had a greater prevalence of HER2-positive cancers than Caucasian women did , . 

We performed an integrated analysis of molecular data from patients through WES, WTS, and mass spectrometry-based proteomics, which provided robust quantification of proteins and phosphoproteins (Fig. ). Whole-exome sequencing of 126 samples revealed 30,193 somatic mutations in the study cohort. Each patient harbored a median of 119 somatic mutations per sample (range: 2–2858 somatic mutations). Significance analysis of the somatic mutations, employing covariates related to the mutational processes, identified the key mutated genes, including TP53, PIK3CA, GATA3, MAP2K4 , and ARID1A (Fig. ). TP53 (60%), PIK3CA (43%), and GATA3 (36%) were the most frequently mutated genes in this cohort. GATA3 somatic mutations were mostly mutually exclusive with TP53 mutations (Fisher’s exact test, P < 0.01), with PIK3CA mutations exhibiting similar mutual exclusivity patterns ( P < 0.01). Overall, approximately 68% of the patients possessed at least one mutation in one of these three genes. The prevalence of somatic mutations in MAP2K4 and ARID1A in breast cancer has also been documented in other studies – , indicating the presence of a congruent mutational profile within the early-onset Asian breast cancer demographic. 

In addition to classical pathological and clinical laboratory-based stratification of the molecular subtypes of the cohorts using hormone receptor and HER2 immunohistochemistry (IHC) analyses, we attempted a robust molecular classification of YBC cohorts based on a multiomics dataset. Unsupervised clustering of the three data platforms (proteome, transcriptome, and phosphoproteome) resulted in discordant stratification unique to each molecular layer. A consensus integrative classification encompassing the three molecular data layers produced five iClusters within the young Korean breast cancer cohort (Fig. ). Cluster 1 was significantly concordant with the classical PAM50-based HER2-enriched subtype. This cluster was characterized by a high prevalence of HER2 amplification and TP53 mutations (69.6% and 82.6% in cluster 1, respectively). The key pathways in this cluster included protein secretion, DNA replication, and repair. Cluster 2 was almost invariably consistent with the basal-like cancer type and presented a high frequency of TP53 mutations and an absence of hormone receptor and HER2 mutations. Cluster 2 was associated with the activation of pathways such as KRAS signaling, autophagy, and cell cycle processes and the negative regulation of apoptosis. Clusters 1 and 2 demonstrated the worst prognosis, with recurrence and metastasis events frequently observed within 5 years after surgical resection (Fig. ). Cluster 3 largely overlapped with the classical luminal B subtype, whereas cluster 4 demonstrated closer alignment with the classical luminal A subtype. The important pathways of cluster 3 and cluster 4 included estrogen signaling, the cell cycle, DNA replication and repair, and fatty acid metabolism (Fig. ). Those clusters of young Korean patients with breast cancer experienced more frequent recurrence and/or death events 5 or more years post-diagnosis of breast cancer (Fig. ). All cases of late relapse in these clusters were ipsilateral (11/11, P < 0.00001, Fisher’s exact test); thus, the late recurrent tumors seemed to have originated from residual and/or resistant tumor cells left from their primary tumors after a sustained period of dormancy rather than from the manifestation of a second primary tumor. In contrast, basal-type and HER2-type Korean YBC demonstrated completely opposite biphasic patterns of survival probability; these groups’ relapse and consequent death events were more clustered in the first 5 years, followed by a plateau in survival probability. The dichotomy of biphasic survival probability between luminal and nonluminal Korean patients with breast cancer is substantially distinct from the clinical course of Caucasian breast cancer, where luminal A breast cancer shows the most indolent phenotype, followed by luminal B breast cancer and HER2 breast cancer, with the basal group having the highest risk patterns . The results of the integrated analysis revealed a distinct subtype, cluster 5, which had a very good prognosis, unlike the other four clusters (Fig. ). The significant pathways of cluster 5 included vasculature development, endothelial cell proliferation, and oxidative phosphorylation. To examine the impact of somatic copy number alterations (SCNAs) on the levels of mRNAs and proteins, we conducted univariate correlation analysis to identify features that exhibited statistically significant correlations (Fig. ). This analysis considered changes in ‘cis’—affecting genes within the altered loci—and ‘trans’—impacting genes located elsewhere. We evaluated the genes in this cohort, each with corresponding SCNA, mRNA, and protein data, by calculating the correlation and its statistical significance for all CNA–mRNA and CNA–protein pairings. Notably, 61% of the CNA-mRNA pairs and 30% of the CNA–protein pairs demonstrated significant positive ‘cis’ correlations. We further scrutinized the subset of these pairs that included well-characterized oncogenes and tumor suppressor genes. ‘Trans’ effects appeared as prominent vertical bands in the analysis, with frequency histograms highlighting significant ‘hot spots’ of ‘trans’ correlation. Interestingly, correlations between CNA and proteins were generally weaker than those observed for CNA and mRNAs. The most pronounced ‘trans’ relationships at the protein level were identified in specific CNA regions, including regions of chromosome 5 and chromosome 8, in early-onset Korean breast cancer. 

Next, we examined the functional impact of genetic alterations on the abundance of both cis -acting (cognate gene products) and trans- acting (other gene products) proteins and their phosphorylation sites (Fig. ). Several cis and trans associations were discernible. These include enhanced expression of HER2 proteins (Fig. ) and strong cis -phosphorylation (Fig. ) in HER2-positive breast cancers, cementing their role as potent drivers of alterations in a subset of breast cancers. Migration and invasion enhancer I ( MIEN1 ) and growth factor receptor-bound protein 7 ( GRB7 ) are genes located proximal to the ERBB2 gene and are concomitantly overexpressed in HER2-positive breast cancer subgroups (Fig. ). MIEN1 has been shown to functionally increase the invasive and migratory phenotypes of various types of cells in breast cancer . Additionally, GRB7 overexpression facilitates the phosphorylation of both AKT and HER2/neu in HER2/neu-overexpressing cells, promoting tumor growth in vivo . These findings suggest additional roles for these proximal genes in breast cancer. Several phosphorylation sites associated with common breast cancer oncogenes were identified and clustered in the cell cycle and DNA repair activities, which is consistent with their established roles in cancer. We investigated which kinase activities were enriched in the molecular subgroups defined in this study (Fig. ). As expected, HER2 kinase activity was upregulated in cluster 1, which was highly correlated with the classical PAM50 Her2 subgroup. Cluster 2 (basal subgroup) was enriched with cell cycle kinases, including CDC7, CDK1, and CDK2. Hormone receptor-positive cancer groups (clusters 3-4) presented increased activity of kinases involved in immune regulation. To identify potential therapeutic targets specific to each iCluster subtype, we used phosphoproteomic data to determine kinase activation (Fig. ). We identified phosphorylated kinases that were enriched in each iCluster subtype via outlier enrichment analysis. The dataset revealed several druggable targets and kinases, including ERBB2 (trastuzumab), in the cluster 1 subtype. The cluster 2 subtype includes isolated PML (tretinoin) and RCSD1 (imatinib). PDCD4 (paclitaxel) was found in cluster 4, whereas CD34 (prednisolone) and TJP1 (dexamethasone) were present in cluster 5. These proteins and kinases are potential targets for subtype-specific treatment. We also investigated the associations between specific genomic alterations and kinase activity (Fig. ). We identified 16 kinases significantly associated with genomic alterations. RB1, PTEN , and MAP2K4 mutations formed a strong subcluster with increased activity of several canonical cell cycle kinases, including the CDK1, CDK2, CDK7, and CDC7 proteins. These associations were also detected in ARID1A- mutant breast cancers. The PLK4, CSNK2A2, HER2, PRKCT, and PRKCH proteins are associated with HER2 genomic alterations. GATA3 genomic alteration was associated with increased activity of the HCK and TYK2 kinases, whereas PIK3CA genomic alteration was linked with the CDK19 kinase. 

We assessed the HER2 driver status by integrating genomic and proteomic data to refine its classification and clinical significance. HER2 status was also grouped by classical IHC and in situ hybridization testing following ASCO-CAP guidelines (HER2+: IHC score of 3+ or an IHC score of 2+ fluorescence in situ hybridization (FISH); silver-enhanced in situ hybridization (SISH): HER2 equivocal status: IHC score of 2+ without FISH results or amplification by FISH without IHC results) . Proteogenomic analysis revealed that 19 patients who presented high levels of both HER2 protein and phosphorylated HER2 were designated as HER2 proteogenomic status positive (HER2 PG+) (Fig. ). HER2 PG+ cases were enriched mostly in cluster 1, followed by cluster 3. The majority of HER2 PG+ samples (16 out of 19) demonstrated largely consistent pathological designations with confirmed HER2-positive expression according to IHC, FISH, and SISH tests. The remaining three patients who were grouped according to HER2 equivocal status demonstrated a high HER2 copy number as a result of droplet digital PCR analysis. Although we identified 28 patients with HER2 gene amplification on the basis of copy number alterations, only 19 of these patients were positive for HER2 by PGstatus. Consistent with the proximal location of the HER2 gene on chromosome 17q12, the copy number, RNA, and protein expression of GRB7 and STARD3 were concomitantly increased in patients with HER2 amplification (Fig. ). The protein levels and phosphorylation metrics of HER2 more robustly and clearly segregated the designation between patients with positive (HER2-PG+) and negative (HER2-PG-) HER2 PG status (Fig. ). In the stratification of patients with breast cancer with HER2 driver dependency, there was a notable incompatibility between conventional PAM50-based classification, classical pathological HER2 type designation, and proteogenomic-based integrative clustering. To gain deeper insights into the biological characteristics that led to sample clustering within the HER2-enriched subtype (HER2E) group despite inconsistent ERBB2 status, we conducted an analysis of phosphosites that are known to canonically relay kinase signaling from the human ERBB2 protein. Patients assigned to the HER2E group were divided into HER2EPG- and HER2EPG+ groups on the basis of their PG status (Fig. ). As anticipated, all PAM50 HER2EPG+ samples presented elevated levels of ERBB2 phosphopeptides. Five additional cases of the PAM50 HER2-enriched subtype were classified into the HER2EPG- subgroup. A close examination of the five cases revealed that these samples shared common activated receptor tyrosine kinase (RTK) downstream pathway signatures, including ARAF, MAPK3, and AKT1 (Fig. ). This finding indicates that the PAM50 simplification scheme mistakenly designates some samples with activated RTKs as belonging to the HER2 subtype, regardless of their subtype-specific nature. Indeed, some samples demonstrated an activated phosphorylation profile of IGF2R, further suggesting that some breast cancer samples might be targetable with precise segregation of RTKs with a suitable driver dependency nature . Quantitative phosphoproteomic analysis revealed that HER2 PG+ cases exhibited elevated levels of HER2 serine phosphopeptides, including serine residues 1054 and 1066, in comparison to PG- cases (Fig. ). Given the close association between the phosphorylation status at serine 1054 (S-1054) and HER2 PG status, we investigated its potential as a predictive biomarker for the sensitivity of breast cancer cell lines to HER2 inhibitors. Despite the limited number of cell lines analyzed, the phosphorylation of S-1054 was found to be correlated with the response to lapatinib and neratinib treatments (Fig. ). We then compared the effects of a HER2 inhibitor on the phosphorylation status of HER2-positive breast cancer cell lines. The phosphorylation level of HER2 at tyrosine 1248 (Y-1248) decreased following treatment with neratinib. The phosphorylation of S-1054 decreased in HER2-positive breast cancer cell lines following neratinib treatment (Supplementary Fig. ). These results suggest that decreased phosphorylation of S-1054 on HER2 could serve as a predictive marker for sensitivity to HER2 inhibitors. Among patients with amplified HER2, the proteogenomic status of HER2 was potentially prognostic, whereas PAM50 HER2 enrichment was not (Fig. ). The patients with HER2 PG+ status tended to have poorer clinical outcomes among the patients with amplified HER2 ( p = 0.024). 

Using transcriptome, proteome, and phosphoproteome datasets, we attempted to classify patients according to their immune signals to understand their tumor immune microenvironments at the molecular and cellular levels (Fig. ). Clusters 1 and 2 were enriched in immune signatures, and the mRNA and protein expression levels of immune checkpoint targets, such as CD276, TIGIT, LAG3, and CTLA4, were also elevated relative to those in the other clusters. Cluster 3 included an increased proportion of patients with fewer immunogenic features, with decreased numbers of T, B, NK, and monocytic lineages, as well as immune checkpoint targets. Multiplex IHC staining of five different immune cell markers (CD8, CD68, FOXP3, PD1, and PD-L1) confirmed the validity of inferring tumor immune microenvironment features from a multiomics dataset (Fig. ). Correlation analyses of multiplexed IHC staining data revealed that YBC immune cells have highly intercorrelated molecular features, confirming the established correlative nature of immune checkpoint molecules, such as augmented PD-L1 staining, which is highly correlated with an increased number of exhausted CD8+ Tregs and immunosuppressive macrophages in the tumor parenchyma (Supplementary Fig. ). CD8-positive T-cell IHC further validated the inferred active immunogenic tumor microenvironment (Fig. ). We also identified three immune-based subtypes with distinct immune characteristics (Fig. ). Among the three clusters, immune-cluster 2 was assigned an immune hot phenotype, immune-cluster 3 an immune cold phenotype, and immune-cluster 1 an immune -intermediate phenotype. There was a positive correlation between the stromal and immune scores according to the three immune-based subtypes (R 2 = 0.4246, p < 0.0001; Fig. ). Classical RNA-seq-based PAM50 classification demonstrated that each subtype of YBC comprises a subset of patients with an active immune signature (Supplementary Fig. ), albeit with higher immune scores in the basal PAM50 subtype (Supplementary Fig. ) (ANOVA, P = 0.007). The active immune tumor microenvironment consisted of elevated levels of stimulatory immune proteins and HLA proteins and increased proportions of tumor-infiltrating lymphocytes (TILs) and myeloid dendritic cells (Supplementary Fig. ). In contrast to the PAM50 stratification, integrative multiomics-based clustering generated immune microenvironment features that segregated cluster 3 from the other clusters (ANOVA, P = 0.001) (Supplementary Fig. ). The evaluation of the clinical characteristics of the poorly immunogenic subgroup tended to be associated with a poor clinical prognosis (Supplementary Fig. ). The clinical association of poor immunogenic features with a poor prognosis was observed in the luminal B PAM50 subtypes in this study cohort (Fig. ), suggesting that immune microenvironment features in specific subgroups of early-onset breast cancer identified through multiomics integration may serve as prognostic molecular markers for stratifying high-risk late-onset luminal B subtype patients. 

To investigate the mutational processes involved in early-onset East Asian breast cancer, we deconvoluted the mutational signatures in the cohort. Four prevalent mutational signatures—BRCAness, APOBEC, mismatch repair defects, and aging signatures—were observed (Fig. ). We also examined germline pathogenic variations, defined as mutations that were reported as the most likely disease-causing variants in ClinVar or mutations that are truncating proteins in 10 genes reported to increase breast cancer susceptibility with a relatively strong penetrance , , . Among the 10 breast cancer susceptibility genes, 10.3% (13/126) of patients harbored pathogenic germline variations, whereas one patient was observed to possess pathogenic variations in both RAD50 and BRCA1. Pathogenic BRCA1/2 germline mutations were more prevalent in Korean patients with YBC, affecting 8.7% of patients in the study cohort compared with 3.8% of patients in the TCGA database (Fisher’s test, P < 0.02). BRCA germline mutations, which are classical biomarkers of BRCAness and PARP inhibitors, only partially explained the proportion of patients with elevated genomic HRD levels (Fig. ). We interrogated the proteome data to search for proteins closely associated with the features of HRD (Fig. ). Among these, 20 highly ranked proteins were chosen whose combination could segregate genomic HRD-positive and HRD-negative cancers (Fig. ). 

Positive correlations existed between the majority of the quantified proteins (>90%) and their corresponding mRNA transcripts within the present study cohort, with a median of 0.34 (Supplementary Fig. ). These correlations were not contingent upon technical factors such as the average protein precursor area or the number of peptide spectral matches. Structural ribosomal proteins and those involved in translational processes presented the lowest correlations between protein and transcript abundance (Supplementary Fig. ), highlighting their function and regulatory control, primarily at the protein level. Notably, proteins derived from transcripts identified in breast cancer drivers, including estrogen receptor responses, exhibited the strongest mRNA‒protein correlations (Supplementary Fig. ), underscoring a direct connection to the protein phenotype in the major canonical breast cancer driver signaling pathway. In hormone receptor-positive breast cancer patients, the subgroup with a high RNA‒protein correlation was associated with a poorer prognosis than the subgroup with a low RNA‒protein correlation (Supplementary Fig. , p = 0.0049). This is reinforced by the greater correlation of genes that respond to estrogen receptor signaling blockade (Supplementary Fig. ). These findings suggest that protein‒RNA correlations could serve as promising prognostic markers for predicting late recurrence in the luminal subtype, necessitating comprehensive external validation accounting for potential confounders in future investigations. 

Patients with YBC have a poor prognosis not only in Western countries , but also in Korea , particularly those with luminal A and luminal B/HER2-type breast cancer . This phenomenon has been attributed to tamoxifen resistance , . Our proteogenomic analysis stratified luminal breast cancer patients into a YBC group with late recurrence (clusters 3 and 4) and very good survival (cluster 5). This classification system can be used to identify patients requiring long-term endocrine therapy. According to the genomic expression analysis, YBC had significantly lower ERα, ERβ, and PR mRNA expression but higher HER-2 and EGFR expression. In addition, lower ERα and higher EGFR mRNA expression were significant predictors of inferior disease-free survival . In our Korean YBC cohort, TP53, PIK3CA and GATA3 were the most frequently mutated genes, which is consistent with the results of the TCGA cohort. The TCGA cohort provides three age groups (≤45, 46-69, and ≥70), and PIK3CA and TP53 were the most common somatic mutations in all three age groups while GATA3 was the third most common and independent somatic mutation in young patients, affecting 15.2% . Proteogenomic analysis integrating the phosphorylation status of HER2 can be used to identify true HER2+ patients. Traditionally, HER2 has been tested for gene amplification or protein overexpression by IHC and ISH. Through proteomic profiling, not only quantitative expression but also posttranslational modification of HER2 could be used for a more precise classification. We identified the proteogenomic status of HER2 as a potential prognostic marker among patients with amplified HER2, as patients with HER2 PG+ status tended to have poorer clinical outcomes than PG- patients. Especially for Asian patients with breast cancer who show a high proportion of HER2 enrichment, HER2 PG status may help predict the clinical outcomes of HER2-targeted therapies. Furthermore, compared with PG- patients, HER2 PG+ patients presented increased phosphorylation of the HER2 protein. Phosphorylation of the HER2 protein has been extensively studied as an essential regulatory mechanism for protein function, but relatively little is known about the effects of serine residues compared with those of tyrosine residues. Our results suggest that the phosphorylation status of the serine 1054 residue is a good marker for estimating the PG status of HER2. PAM50 is a classification tool used to characterize and predict the outcomes of patients with breast cancer – . However, our results demonstrated that the proteogenomic status of HER2 was a potential prognostic factor, whereas PAM50 HER2 enrichment was not. Patients with a HER2 PG+ status tended to have poorer clinical outcomes. The HER2-positive subtypes show ethnic disparities, and Korean patients with breast cancer have a high degree of HER2 enrichment , , . Our results suggest that HER2 PG status is required to select true-HER2+ patients who respond well to anti-HER2 treatment. As expected, all of the PAM50 HER2EPG+ samples had high levels of HER2 phosphopeptides, whereas the PAM50 HER2EPG- samples had markedly lower levels but elevated levels of phosphorylation of other ERBB family members and of the mitogen-activated protein kinase (MAPK) signaling pathway compared with PG+ samples. These findings suggest that complementary pathways that can replace HER2 signaling may be targetable in PAM50 HER2E tumors without HER2 amplification. The tumor microenvironment (TME) is a critical factor that influences the response to immunotherapy. The cytotoxic subset of CD8+ T cells can eliminate cancer cells and is associated with improved patient survival. Conversely, immunosuppressive regulatory CD4+ T cells (Tregs) or macrophages are associated with a poorer prognosis . “Hot” tumors characterized by high levels of T-cell infiltration generally exhibit more favorable responses to immunotherapies than “cold” tumors, which have limited T-cell infiltration . The existence of TILs within the TME indicates the presence of an inherent immune response against the tumor, which is linked to a more favorable prognosis and increased responsiveness to chemotherapy . Compared with more immunogenic cancers, breast cancer, which is traditionally considered to be poorly immunogenic, or “cold” with a modest tumor mutation burden (TMB) and a low TIL count, shows a notable prevalence of APOBEC mutational signatures – . The mutational load in BC exhibits significant variability, and tumors with higher TMB may exhibit a more favorable response to immune checkpoint inhibitors (ICIs) . Although PD-L1 expression is commonly used as a biomarker, it alone is not sufficient as a predictor, demonstrating its predictive value mainly in metastatic TNBC. TNBC and HER2-positive BC are considered more immunogenic, whereas hormone receptor-positive BC is considered less immunogenic. However, certain patients with HER2-positive or hormone receptor-positive BC may also show promise for immunotherapy in addition to TNBC . The study cohort revealed a clinical association between an unfavorable prognosis and poor immunogenic characteristics, specifically within the luminal B PAM50 subtypes. This implies that, in cases where hot immune groups are identified within specific subtypes, it is essential to consider immunotherapeutic interventions. Our proteome-based analysis identified HRD-associated proteins and an HRDness that could facilitate the adaptive and versatile identification of patients with YBC who may benefit from treatment with PARP inhibitors. Additionally, we demonstrated that protein–RNA correlations can be used to predict the late recurrence of hormone receptor-positive breast cancer. These findings have valuable clinical implications for improving survival outcomes and reducing treatment-related toxicity by guiding the optimal duration and selection of adjuvant treatment for patients with breast cancer. These findings can be clinically applied to guide the optimal duration and selection of adjuvant treatment for patients with breast cancer, which is particularly crucial for patients diagnosed with breast cancer at a young age. Within each standard molecular subtype of breast cancer, we identified functionally significant protein groups whose differential abundance closely correlated with the clinical progression of breast cancer. Moreover, we derived a recurrence-predictive cluster capable of predicting late recurrence, specifically for the luminal subtype, and validated its efficacy in various patient cohorts. Prognostication currently used for breast cancer depends largely on clinical subtypes and mRNA-based information. This proteome-integrated prognostic index plays a crucial role in guiding treatment duration decisions for early-onset breast cancer, thereby contributing to improved patient stratification and personalized treatment approaches. 

Supplementary information 

