A molecular classification system for estimating radiotherapy response and anticancer immunity for individual breast cancer patients

Objective Radiotherapy is a cornerstone of breast cancer therapy, but radiotherapy resistance is a major clinical challenge. Herein, we show a molecular classification approach for estimating individual responses to radiotherapy Methods Consensus clustering was adopted to classify radiotherapy-sensitive and -resistant clusters in the TCGA-BRCA cohort based upon prognostic differentially expressed radiotherapy response-related genes (DERRGs). The stability of the classification was proven in the GSE58812 cohort via NTP method and the reliability was further verified by quantitative RT-PCR analyses of DERRGs. A Riskscore system was generated through Least absolute shrinkage and selection operator (LASSO) analysis, and verified in the GSE58812 and GSE17705. Treatment response and anticancer immunity were evaluated via multiple well-established computational approaches. Results We classified breast cancer patients as radiotherapy-sensitive and -resistant clusters, namely C1 and C2, also verified by quantitative RT-PCR analyses of DERRGs. Two clusters presented heterogeneous clinical traits, with poorer prognosis, older age, more advanced T, and more dead status in the C2. The C1 tumors had higher activity of reactive oxygen species and response to X-ray, proving better radiotherapeutic response. Stronger anticancer immunity was found in the C1 tumors that had rich immune cell infiltration, similar expression profiling to patients who responded to anti-PD-1, and activated immunogenic cell death and ferroptosis. The Riskscore was proposed for improving patient prognosis. High Riskscore samples had lower radiotherapeutic response and stronger DNA damage repair as well as poor anticancer immunity, while low Riskscore samples were more sensitive to docetaxel, doxorubicin, and paclitaxel. Conclusion Our findings propose a novel radiotherapy response classification system based upon molecular profiles for estimating radiosensitivity for individual breast cancer patients, and elucidate a methodological advancement for synergy of radiotherapy with ICB.


Introduction
Breast cancer is a significant global health concern, accounting for a substantial number of cancer-related deaths among women worldwide.The latest global cancer burden data released by the World Health Organization's International Agency for Research on Cancer (IARC) for 2020 shows that there are 2.26 million new breast cancer cases worldwide (1).About 13% of women are likely to be diagnosed with breast cancer in their lifetime, and the rate of increase is 0.3% per year, which seriously endangers the health of patients (2).
Despite advancements in early detection and treatment strategies, a considerable proportion of breast cancer patients experience disease recurrence and metastasis, leading to poor clinical outcomes.As a fundamental component of breast cancer therapy, 75% of breast cancer patients will receive radiation therapy, which is employed to eradicate subclinical tumor lesions and reduce the risk of local recurrence (3).Research data suggest that radiotherapy can improve the local control of breast cancer by 60% to 70%, and increase the absolute survival rate by 10% (3).However, radiotherapy resistance remains a major clinical challenge, limiting its efficacy and compromising patient outcomes (4).
In recent years, molecular profiling techniques have revolutionized cancer research by providing a comprehensive view of the molecular landscape of tumors (5).According to the status of estrogen receptor (ER), progesterone receptor (PR), ki-67, and HER-2 expression, breast cancer can be classified into four molecular subtypes: Luminal A, Luminal B, HER-2 over-expression, and triple-negative breast cancer.This has become a model for the precision treatment of breast cancer (6).Although this molecular subtyping has been extensively studied in the context of prognosis and treatment selection, there is a paucity of research focusing specifically on molecular classification for estimating radiotherapy response in early breast cancer patients (7).
Insight of the molecular mechanisms underlying radiotherapy response in breast cancer is crucial for improving treatment outcomes and personalizing patient management.Traditionally, clinical factors such as tumor stage, histological grade, and hormone receptor status have been used to guide treatment decisions.However, these factors often fail to accurately predict individual responses to radiotherapy.With the progress of molecular medicine radiotherapy technology, the radiotherapy of breast cancer is also developing towards the direction of individualization.Revealing the radioresistance of breast cancer by molecular mechanism has become a hot topic in radiotherapy research.Radioresistance is a complex process that is generally associated with radiation-induced DNA damage repair, cell cycle dysregulation, cancer stem cell properties, and epithelialmesenchymal transition (8).Therefore, there is an increasing interest in identifying molecular markers and developing molecular classification systems that can better evaluate radiotherapy response in breast cancer patients.The ability to accurately predict individual responses to radiotherapy would help identify patients who are likely to benefit from this treatment modality and exclude others from potential side effects of radiotherapy.
In recent years, immunotherapy has become a hot research topic for scholars at home and abroad as a new treatment.Clinical trials have shown that PD-1/PD-L1 inhibitors can improve the antitumor efficacy in triple-negative breast cancer (9).Radiotherapy may lead to the overexpression of PD-L1 on tumor cells by activating PI3K/AKT, signal transduction, and transcription factor activation.Moreover, PD-L1 stimulates cell migration and promotes the process of epithelial-mesenchymal transition, thereby inducing radioresistance (10).In the basic research of breast cancer, it has been found that PD-1/PD-L1 inhibitors combined with radiotherapy have a stronger anti-tumor effect than single-mode treatment, and the survival time of mice is significantly higher than that of the control group (11,12).
The objectives of this study are to classify breast cancer patients into radiotherapy-sensitive and -resistant clusters and assess the clinical and molecular characteristics of these clusters.In this study, we employed consensus clustering analysis to identify distinct clusters of patients based on the expression profiles of radiotherapy response-related genes.The stability of the classification was validated in independent patient cohorts.Furthermore, we developed a Riskscore system using LASSO analysis to improve patient prognosis.Additionally, we evaluated treatment response and anticancer immunity using well-established computational approaches.Our findings provide a methodological advancement for the synergy of radiotherapy with immune checkpoint blockade (ICB), potentially enhancing treatment outcomes for breast cancer patients.

Data acquisition and processing
A total of 3693 radiotherapy response-related genes (RRGs) were acquired from previously published literature (13).Supplementary Table 1 summarizes the RRG list.RNA sequencing data of The Cancer Genome Atlas (TCGA) Breast Cancer (TCGA-BRCA) were gathered from the Genomic Data Commons (https:// portal.gdc.cancer.gov/)utilizing the TCGAbiolinks package (14).Following the removal of samples without complete prognostic information, 995 samples were included.The raw read count data were transformed to transcripts per kilobase million, with further log-2 conversion.From the Gene Expression Omnibus (https:// www.ncbi.nlm.nih.gov/geo/),two independent breast cancer cohorts: GSE58812 (n=107) (15) and GSE17705 (n=298) (16) were utilized for external verification.The two microarray cohorts were based on the Affymetrix platform, and background correction and standardization were implemented via a robust multiarray averaging approach utilizing the affy package (17).Basic information and demographic data of TCGA-BRCA, GSE58812, and GSE17705 datasets were summarized in Supplementary Table 2.

Selection of differentially expressed RRGs
Differential expression analysis on RRGs was conducted by comparing breast cancer patients who received radiotherapy and those who did not receive radiotherapy.This analysis was achieved via limma package (18).DERRGs were identified with adjusted p-value<0.05.

Genomic variation analysis
Copy number variation (CNV) data were gathered from the Fire Browse (http://firebrowse.org/).The GISTIC2.0 computational method was adopted for the estimation of gene gains and losses (19).Somatic mutation profiling was obtained from cBioPortal (https://www.cbioportal.org/)and was analyzed based on the matfools package (20).

Functional enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were achieved utilizing clusterProfiler (21).The gene sets of reactive oxygen species (ROS), response to X-ray and DNA damage repair pathways were gathered from the Molecular Signatures Database (22).The enrichment score was quantified via single-sample gene set enrichment analysis (ssGSEA) (23).

Consensus clustering analysis
Univariate-cox regression analysis on DERRGs with prognosis was carried out.DERRGs with p-value<0.05were employed for consensus clustering analysis via the ConsensusClusterPlus package (24).The optimal number of clusters and stability were determined based on the consensus matrix heatmap, the proportion of ambiguous clustering (PAC), and principal component analysis (PCA).The repeatability and accuracy of the classification were validated via nearest template prediction (NTP) analysis (25).This analysis was conducted based on the top 100 upregulated markers of each cluster in the GSE58812 cohort using the CMScaller package (26).

Least absolute shrinkage and selection operator analysis
Using the glmnet package (38), LASSO analysis was implemented based on prognostic DERRGs.After identifying the minimum lambda, DERRGs with coefficient ≠0 were further selected through predict.cv.glmnet function.Riskscore was subsequently defined through a combination of expressions of the selected DERRGs with given coefficients.Patients were stratified into low-and high-Riskscore groups with the median Riskscore.The reliability and stability of the model were proven in the GSE58812 and GSE17705 datasets.

Nomogram establishment
Uni-and multivariate-cox regression analyses on RiskScore, and clinical variables with TCGA-BRCA prognosis were implemented.Based upon independent prognostic variables, a nomogram was generated utilizing the rms package.Additionally, calibration curves were plotted for comparing the nomogramestimated and actual outcomes.

Prediction of transcription factors
By adopting the NetworkAnalyst online tool (http:// www.networkanalyst.ca)(39), a transcription factor-gene interaction network was visualized.

Drug sensitivity analysis
Half-maximal inhibitory concentration (IC50) of chemotherapy agents was computed via pRRophetic analysis (40) based upon the Genomics of Drug Sensitivity data (41).

Drug-gene interaction analysis
Drugs that potentially targeted the DERRGs from the Riskscore were estimated based on the Drug-Gene Interaction Database (www.dgidb.org)(42), and a drug-gene interaction network was generated utilizing the Cytoscape tool (43).
The breast cancer samples were classified into the resistant and sensitive groups based on their response to radiotherapy, where patients with complete response (CR), partial response (PR), and stable disease (SD) comprised the radiotherapy sensitive group (Rad-S, n=10) and progressive disease (PD) belonged to the radiotherapy resistant group (Rad-R, n=8).The clinical pathological characteristics of the patients are listed in Supplementary Table 3.
The tissues were obtained from breast cancer patients participating in a previous clinical study in Tongji Hospital.The ethics approval was not required by the Ethics Committee of Tongji Hospital Affiliated with Tongji Medical College of Huazhong University of Science and Technology, because the clinical study was previously approved by the Ethics Committee.

Statistical analysis
Data were analyzed via appropriate R packages (version 3.6.1).Continuous variables between two groups were analyzed with Student's t-test or Wilcoxon test, while categorical data were evaluated via the chi-square test.Pearson's test was used for correlation analysis.Survival analysis was executed through Kaplan-Meier curves and log-rank test utilizing the survival package.The specificity and sensitivity were appraised via receiver operating characteristic curves (ROCs), and the area under the curve (AUC) was quantified utilizing the pROC package.P-value<0.05was statistically significant.

Classification of breast cancer patients into radiotherapy-sensitive and -resistant consensus clusters
TCGA-BRCA samples were classified as two consensus clusters based on the transcriptome values of the prognostic DERRGs (Figure 2A).PAC score was relatively small at cluster number k=2 (Figure 2B).This demonstrated that the optimal k value was 2. PCA also unveiled the diverse transcriptome profiling in two clusters (Figure 2C).Many radiotherapy-sensitive genes were up-regulated in the C1 tumors, while radiotherapy-resistant genes were upregulated in another cluster (Figure 2D).Thus, we inferred the C1 tumors were regarded as radiotherapy-sensitive clusters, while the C2 tumors were regarded as radiotherapy-resistant clusters.Worse survival was detected in the C2 versus C1 patients (Figure 2E).The GSE58812 dataset was adopted to prove the classification.The top 100 up-regulated marker genes in each consensus cluster were selected (Figure 2F; Supplementary Table 5), and samples with p-value<0.05were extracted for quantification and assessment of prediction confidence (Figure 2G).Consistently, the C2 posessed a less favorable prognosis than the C1 (Figure 2H).In addition, the mRNA levels of CIITA and IL27RA are significantly downregulated, while the mRNA levels of ZFP41 and N4BP3 are significantly increased in radiotherapy-sensitive tumors, rather than radiotherapy-resistant tumors (Figure 2I).The detection results are consistent with the predictions of this system.Hence, the classification was reliable and repeatable.

Two consensus clusters are characterized by different clinical traits and radiotherapy responses
The C2 had more cases with age ≥65, more advanced T, and dead status in comparison to the C1, but without difference in N, M, and stage (Figures 3A-F).Radiotherapy results in DNA damage directly through ionization or indirectly through generating ROS, thus destroying tumor cells (44).The C1 tumors posessed stronger ROS activity and responses to X-ray (Figures 3G, H), further demonstrating that patients in this cluster better responded to radiotherapy.

Two consensus clusters present heterogeneous anticancer immunity and immune escape
Most immune cell populations (e.g., natural killer cells, natural killer T cells, activated dendritic cells, activated CD4 and CD8 T cells) displayed more infiltration in the C1 versus another cluster, demonstrating stronger anticancer immunity in C1 tumors (Figure 5A).Lower exclusion, dysfunction, IFNG, and TIDE scores were found in C2 tumors, inferring patients in the cluster potentially responded to ICB (Figures 5B-E).In addition, the C1 exhibited a similar expression profile to that of samples responding to the PD-1 antibody (Figure 5F), further proving stronger anticancer immunity in C1 tumors.Immunogenic cell death that can be induced by radiotherapy exerts a crucial role in evoking systemic immune response against tumors (45).Many immunogenic cell death molecules, e.g., TLR2/3/4/7/9, CALR, CGAS, CLEC4E, CLEC7A, DDX58, FPR1/2, HMGB1, IL33, NLPR3, and P2RX7 were notably up-regulated the C1 (Figure 5G), revealing activated immunogenic cell death in the cluster.Ferroptosis has been implicated in anticancer immunity and immune response (46).The C1 presented remarkable up-regulation of many ferroptosis molecules, such as LSP1, CAT, SOD2, PRNP, FTL, FES, GLRX, PFKP, PDIM1, MBP, HHEX, and GPX3 (Figure 5H).Overall, the C1 tumors were characterized by potent anticancer immunity.6A).Using the formula, the RiskScore of each TCGA-BRCA sample was computed.Under the median RiskScore, TCGA-BRCA samples were stratified into two groups (Figure 6B).Survival difference was found in two groups, with poorer survival probability for the high RiskScore group versus another group (Figure 6C).ROCs were plotted to appraise the prediction efficiency.The one-, three-, and five-year AUC values were 0.785, 0.800, and 0.800, respectively, proving that the RiskScore was capable of estimating patient prognosis.The excellent efficacy was externally demonstrated in the GSE17705 and GSE58812 datasets (Figures 6D, E).
Uni-in combination with multivariate-cox regression results demonstrated the independency of the RiskScore in prognostication in addition to age and stage (Figures 6F, G).Based upon them, a nomogram was generated for estimation of one-, three-, five-, eight-, and ten-year survival probability (Figure 6H).To verify its performance, calibration curves were conducted.Consequently, the nomogram-estimated prognosis was close to the actual outcomes (Figure 6I).Overall, our nomogram provided a method to estimate survival outcomes via integration of the RiskScore, age, and stage.

The RiskScore positively correlates to radiotherapy resistance
The RiskScore was negatively connected to response to X-ray, especially for the C2 tumors (Figure 7B).A positive relationship between the RiskScore and ROS was also found in the C2 tumors (Figure 7C).In addition, the RiskScore was positively linked with DNA damage response pathways, e.g., base excision repair, direct reversal of damage, homologous recombination, and nucleotide excision repair (Figure 7D).Altogether, the RiskScore was positively connected to radiotherapy resistance.

The RiskScore is negatively associated with anticancer immunity
The study executed multiple well-established computational methods to estimate the abundance of immune cells.Intriguingly, most immune cells, notably cytotoxic Tlymphocytes that are critical effectors of anticancer immunity (47), presented higher infiltration in low-RiskScore tumors (Figure 8A).Moreover, the RiskScore was negatively linked with cytotoxic T-lymphocytes (Figure 8B), indicative of a negative association of the RiskScore with anticancer immunity.We also focused on the remarkably negative relationships of the RiskScore with common immune checkpoint molecules, e.g., PDCD1 (PD-1) and its ligands CD274 (PD-L1) and PDCD1LG2 (PD-L2) (Figure 8C) (48).This indicated that low-RiskScore tumors displayed a stronger sensitivity to ICB.The RiskScore was also found to be negatively connected to most immunogenic cell death molecules, such as IL33, TLR2/3/4/7, CLEC4E, NLRP3, CGAS, and FPR1/2 (Figure 8D), demonstrating a negative association of the RiskScore with immunogenic cell death.Docetaxel, doxorubicin, and paclitaxel are routinely applied in chemotherapy for breast cancer.Here, we evaluated the difference in sensitivity to these chemotherapeutic agents.Consequently, the low RiskScore group presented significantly lower IC50 values of docetaxel, doxorubicin, and paclitaxel versus the high RiskScore group (Figures 9A-C).Hence, low-RiskScore patients were more sensitive to the above chemotherapeutic drugs.
Small molecules that potentially targeted the DERRGs from the RiskScore model were also analyzed.As illustrated in Figure 9D,

Discussion
The present study has provided a comprehensive molecular classification system for predicting radiotherapy response and anticancer immunity in individual breast cancer patients.This classification system, based on the consensus clustering of differentially expressed radiotherapy response-related genes (DERRGs), has identified two distinct clusters, C1 and C2, with C1 tumors demonstrating a higher sensitivity to radiotherapy and stronger anticancer immunity.The identification of these two distinct clusters has significant implications for the personalized treatment of breast cancer patients, as it helps the prediction of individual responses to radiotherapy, thereby enabling the decision of treatment strategies to enhance therapeutic efficacy and minimize adverse effects.
The C1 cluster, characterized by a higher activity of Reactive Oxygen Species (ROS) and a stronger response to X-ray, demonstrated a better radiotherapeutic response.ROS accumulation can lead to DNA damage and genomic instability, induce multiple forms of cell death, and play an important role in cancer development and progression (49).In preclinical studies of breast cancer, ROS-induced nanophototherapy systems have been shown to increase radiosensitivity (50).This finding is consistent with previous studies (51) that have highlighted the role of reactive oxygen species in mediating the cytotoxic effects of radiotherapy.Furthermore, the C1 tumors exhibited a stronger anticancer immunity, as evidenced by the rich immune cell infiltration and the similar expression profiling to patients who responded to anti-PD-1.In this study, we validated the high expression of MHC class II transactivator (CIITA) in the C1 population through quantitative PCR.CIITA is a master controller of antigen presentation and a crucial factor in adaptive immunity.In recent years, CIITA has been speculated to be an important exploration target in anti-checkpoint blockade immunotherapy and anti-tumor vaccination (52,53).Therefore, the C1 tumors may be more sensitive to ICB, which has emerged as a promising therapeutic strategy for various types of cancer.However, the effects of RGGs on chemotherapy response in the C1 population are somewhat contradictory, with some genes enhancing chemotherapy sensitivity (IL7R) (54) while others promoting chemotherapy drug resistance (PLAC8, CCL5) (55,56).
In contrast, the C2 cluster was associated with poorer prognosis, older age, more advanced T stage, and more dead status.These tumors demonstrated a lower response to radiotherapy and a stronger DNA damage repair capacity, suggesting a higher resistance to radiotherapy.Reactive oxygen species (ROS) are one of the major inducers of DNA damage, and one of the key determinants of the efficacy of anti-tumor therapy is the severity of DNA damage resulting in tumor cells.However, cancer cells are able to adapt DNA repair pathways to make tumors resistant to treatment (57).Inhibitors targeting the DNA repair system may enhance the radiosensitivity of tumors, but a few studies have suggested that the loss of DNA repair function may increase the rate of gene mutation and eventually lead to radioresistance.Therefore, the relationship between DNA damage repair and radiotherapy resistance is still controversial and needs further study (58).This finding underscores the need for alternative therapeutic strategies for patients with C2 tumors, such as the use of DNA damage repair inhibitors to sensitize these tumors to radiotherapy.Moreover, the C2 tumors exhibited poor anticancer immunity, which highlights the need for further research to elucidate the mechanisms underlying the immune escape.
The RiskScore system proposed in this study represents a significant advancement in the prognostication of breast cancer patients.This system, based on the expression of DERRGs, was capable of predicting patient prognosis with high accuracy, as evidenced by the high area under the curve (AUC) values.Importantly, the RiskScore system was independent of age and stage, suggesting that it could provide additional prognostic information beyond these traditional prognostic factors.Furthermore, the RiskScore system was associated with radiotherapy response and anticancer immunity.The high-RiskScore samples demonstrate lower radiotherapeutic response and poor anticancer immunity, while low-RiskScore samples are more sensitive to chemotherapy drugs such as docetaxel, doxorubicin, and paclitaxel.This suggests the RiskScore system could potentially help the decision of therapeutic strategies for individual patients.
In conclusion, this study has proposed a novel radiotherapy response classification system based on molecular profiles for evaluating radiosensitivity in individual breast cancer patients.This system, combined with the RiskScore system, could potentially revolutionize the personalized treatment of breast cancer by enabling the prediction of individual responses to radiotherapy and immunotherapy.However, further validation in prospective clinical trials is warranted to confirm the clinical significance of these systems.Moreover, future research should focus on elucidating the molecular mechanisms underlying the differential radiosensitivity and anticancer immunity of the identified clusters, which could guide the development of new therapeutic targets.

1
FIGURE 1 Multi-omics analysis of DERRGs in breast cancer.(A, B) Identification of DERRGs through comparing patients with radiotherapy with those without radiotherapy.(C) Somatically mutated DERRGs across breast cancer.Mutation forms are marked by unique colors, while DERRGs are ranked by mutated frequency.(D) Copy number gains and losses of DERR © (E) Univariate-cox regression results on DERRGs with patient survival.(F-I) GO and KEGG pathways enriched by DERRGs.

2 3 4
FIGURE 2 Classification of breast cancer patients into radiotherapy-sensitive and -resistant consensus clusters.(A) Consensus matrix heatmap at cluster number k=2.White means samples never get together; blue means the samples are always clustered together.(B) PAC distribution at cluster numbers k=2~9.(C) PCA verifying the different transcriptome profiles in two clusters.(D) Expression values of the prognostic DERRGs in two clusters.Blue to red denotes down-to up-regulated expre © on.(E) Survival probability of two clusters.(F) The top 100 up-regulated marker genes in each cluster.(G, H) Verification of the patient classification, and survival difference via NTP method in the GSE58812 cohort.(I) Comparison of CIITA, 27RA, ZFP41, and N4BP3 expression in breast cancer tissues that are sensitive (n=10) or resistant (n=8) to radiotherapy by quantitative real-time PCR detection.

5 6
FIGURE 5 Two consensus clusters exhibit heterogeneous anticancer immunity and immune evasion mechanisms.(A) Abundance of immune cell types in two clusters.(B-E) Different exclusion, dysfunction, IFNG, and TIDE scores in two clusters.(F) Submap for investigating the similarity in the expression profiling of two clusters with that of samples that responded or did not respond to PD-1 or CTLA4 antibody.(G, H) Difference in immunogenic cell death and ferroptosis molecules between two clusters.

7
FIGURE 7    Transcriptional regulatory mechanisms underlying the RiskScore and a positive association of the RiskScore with radiotherapy resistance.(A) Transcriptional factor-gene network.Rhombus, transcriptional factor; circle, the DERRG from the RiskScore.(B, C) Relationships of the RiskScore with response to X-ray and ROS in the C1 and C2 tumors.(D) Correlation between the RiskScore with DNA damage repair pathways.The bubble from green to pink denotes the association from negative to positive.P-value is represented by a unique line color.

8
FIGURE 8 Association of the RiskScore with immunogenicity.(A) Distribution of the infiltration of distinct immune cell populations in low-and high-RiskScore groups by use of multiple algorithms.Green to red denotes weak to strong infiltration.(B) Association of the RiskScore with infiltrative immune cells.The bubble location indicates the correlation coefficient.(C) Correlation between the RiskScore with known immune checkpoints.The bigger the bubble, the stronger the correlation.(D) Relationships of the RiskScore with immunogenic cell death molecules.The bubble from green to pink shows the correlation from negative to positive.P-value is represented by a unique line color.

9
FIGURE 9 Analysis of drug sensitivity and druggable targets.(A-C) Estimated IC50 values of docetaxel, doxorubicin, and paclitaxel in low-and high-RiskScore groups.The larger the IC50, the weaker the sensitivity to a drug.(D) Drug-gene network.Triangle, drug; circle, the DERRG from the RiskScore model.