A Multiomics Signature Highlights Alterations Underlying Homologous Recombination Deficiency in Triple-Negative Breast Cancer

Homologous recombination (HR) is a key pathway in DNA double-strand damage repair. HR deficiency (HRD) occurs more commonly in triple-negative breast cancers (TNBCs) than in other breast cancer subtypes. Several clinical trials have demonstrated the value of HRD in stratifying breast cancer patients into distinct groups based on their responses to poly(ADP ribose) polymerase inhibitors and chemotherapy. We retrospectively collected TNBC samples to establish a multiomics cohort (n = 343) and explored the biological and phenotypic mechanisms underlying the better prognosis of patients with high HRD scores. Gene set enrichment analysis was conducted to elucidate the underlying pathways in patients with low HRD scores, and a radiomics model was established to predict the HRD score via a noninvasive method. Multivariable Cox analysis revealed the independent prognostic value of a low HRD score (hazard ratio 2.20, 95% confidence interval 1.05–4.59; p = 0.04). Furthermore, amino acid and lipid metabolism pathways were highly enriched in tumors from patients with low HRD scores, which was also demonstrated by differential abundant metabolite analysis. A noninvasive radiomics method was developed to predict the HRD status and it performed well in the independent validation cohort (support vector machine model: area under the curve [AUC] 0.739, sensitivity 0.571, and specificity 0.824; logistic regression model: AUC 0.695, sensitivity 0.571, and specificity 0.882). We revealed the prognostic value of the HRD score, predicted the HRD status with noninvasive radiomics features, and preliminarily explored druggable targets for TNBC patients with low HRD scores.

HRD in stratifying breast cancer patients into distinct groups based on their responses to poly(ADP ribose) polymerase inhibitors and chemotherapy. Methods. We retrospectively collected TNBC samples to establish a multiomics cohort (n = 343) and explored the biological and phenotypic mechanisms underlying the better prognosis of patients with high HRD scores. Gene set enrichment analysis was conducted to elucidate the underlying pathways in patients with low HRD scores, and a radiomics model was established to predict the HRD score via a noninvasive method. Results. Multivariable Cox analysis revealed the independent prognostic value of a low HRD score (hazard ratio 2.20, 95% confidence interval 1.05-4.59; p = 0.04). Furthermore, amino acid and lipid metabolism pathways were highly enriched in tumors from patients with low HRD scores, which was also demonstrated by differential abundant metabolite analysis. A noninvasive radiomics method was developed to predict the HRD status and it performed well in the independent validation cohort (support vector machine model: area under the curve [AUC] 0.739, sensitivity 0.571, and specificity 0.824; logistic regression model: AUC 0.695, sensitivity 0.571, and specificity 0.882). Conclusions. We revealed the prognostic value of the HRD score, predicted the HRD status with noninvasive radiomics features, and preliminarily explored druggable targets for TNBC patients with low HRD scores.
Triple-negative breast cancer (TNBC) is a subtype of breast cancer defined by the absence of the estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2). 1,2 TNBC typically occurs in young patients and progresses rapidly, leading to poor prognosis. 3,4 The lack of effective therapeutic targets indicates an urgent need to identify novel targets and screen patients who are suitable for specific therapies. In recent years, the development of sequencing techniques has highlighted the heterogeneity of breast cancer subtypes with distinct biological characteristics. [5][6][7] Our previous study classified TNBC into four subtypes, luminal androgen receptor (LAR), basal-like immune-suppressed (BLIS), immunomodulatory (IM), and mesenchymal-like (MES), providing practicable treatment strategies. [5][6][7] DNA in cellular nuclei is continuously exposed to damage risk factors, causing a variety of patterns of lesions, including single-strand breaks (SSBs), doublestrand breaks (DSBs), base mismatches, insertions, and deletions. 8 Among them, DSBs are typically repaired by homologous recombination (HR) and nonhomologous end joining (NHEJ). 9 Cancer cells harboring mutations in a gene orchestrating HR, such as the BRCA1, BRCA2, RAD51, and PALB2 genes, are more vulnerable to DNAdamaging drugs such as platinum, and to the SSB repair suppressor poly(ADP ribose) polymerase (PARP) inhibitors. [10][11][12][13] Notably, Timms et al. introduced the HR deficiency (HRD) score in a retrospective study by summing the telomeric allelic imbalance (TAI), loss of heterozygosity (LOH), and large-scale state transition (LST) values. 14-17 Telli et al. defined HRD based on a score C 42 or BRCA 1/2 mutation and found that breast cancer patients with HRD were prone to achieving a pathological complete response (pCR) to neoadjuvant chemotherapy. In addition, a high HRD score (HRD C42) but BRCA 1/2 wild-type was positively correlated with a higher pCR rate, further demonstrating the value of the HRD score in breast cancer. 18 Data from The Cancer Genome Atlas (TCGA) revealed that germline BRCA1/2 mutations exist in approximately 10% of patients with TNBC; [19][20][21] however, the mechanism underlying the better clinical outcomes of high HRD TNBC patients has not been thoroughly investigated.
In this study, we integrated multiomics information from our TNBC cohort at Fudan University Shanghai Cancer Center (FUSCC). Genomics, transcriptomics, metabolomics, and radiomics were utilized to calculate the HRD scores, analyze the differences between patients with high and low HRD scores, and noninvasively predict the HRD status.

Study Population
A total of 343 TNBC patients with preoperative magnetic resonance imaging (MRI), primary tumor tissue, and blood samples obtained from FUSCC made up the final study group. All patients underwent surgery at FUSCC between 2007 and 2014, and the median follow-up duration was 47.3 months. A total of 67 patients experienced relapse and/or metastatic progression during this period. The clinical data included patient age, tumor size, lymph node status and Ki-67 index, adjuvant chemotherapy regimen, and metastatic size.
In this study, all 343 patients were included in the HRD analysis. Among them, 118 patients were available for radiomics analysis and had preoperative breast contrastenhanced (CE) MRI, and 252 sections stained with hematoxylin and eosin (H&E) and immunohistochemistry (IHC) staining data were available. In addition, 252 samples with transcriptomic sequencing data and 246 samples with polar metabolomics and lipidomics data were also eligible for analysis. The generation of transcriptomics and metabolomics data has been described previously. 22,23 All the samples from TNBC patients were used for (1) HRD analysis and prognosis assessment (n = 343); (2) radiomics signature analysis and validation for HRD identification (n = 118); and (3) molecular and biological assessments of associations between HRD and other multiomics, including pathology, transcriptomics, and metabolomics (n = 252, 252, and 246, respectively). Multiomics detection and pathological indices assessment were conducted using surgical biopsy samples as described in a previous study. 23 This study was approved by the Medical Ethics Committee of FUSCC, and written informed consent was provided by each patient. The study was conducted in accordance with the Declaration of Helsinki as revised in 2013.

Homologous Recombination Deficiency (HRD) Analysis
As summarized in a previous study, 14 the HRD score was calculated as the arithmetic sum of three scores (TAI, LOH, and LST) using segments with integer copy numbers produced by ASCAT. The methods used to calculate these three scores have been described in a previous study. 14 In brief, the TAI score was defined as the number of regions with allelic imbalance that were longer than 11 Mb and extended to one of the subtelomeres but did not cross the centromere. The LOH score was defined as the number of LOH regions longer than 15 Mb but shorter than the whole chromosome, whereas the LST score was defined as the number of break points between regions longer than 10 Mb after filtering out regions shorter than 3 Mb. To diminish the effect of ploidy, the LST score was modified according to the formula LSTm = LST-kP, where P is ploidy and k is a constant of 15.5. The HRD cohort was divided into high (HRD score C42) and low (HRD score \42) HRD cohorts. 18

Analysis of Differentially Abundant Metabolites and Differentially Expressed Genes
The differentially expressed genes were determined through the Limma R package, 24 and the differential expression analysis outputs of Limma were used to generate the ranked list file. Gene set enrichment analysis (GSEA) was performed using the clusterProfiler R package, 25 and the differential abundance of metabolites was determined by calculating the Mann-Whitney U values for all the detected metabolites.

Magnetic Resonance Imaging Technique and Radiomics Analysis
MRI scans were obtained using a 1.5-T dedicated spiral breast MRI system (Aurora Imaging Technology, Aurora Systems, Inc., Vancouver, BC, Canada) described elsewhere. In this study, CE-MRI data were used for radiomic analysis. The T1-weighted fat-suppressed sequence (time of echo/time of repetition 4.8/29 ms; thickness 1.1 mm; field of view 360 9 360 mm; matrix 360 9 360) included one pre-contrast image and four post-contrast images following a bolus injection of gadopentetate dimeglumine at 2 mL/s (Magnevist, Bayer Schering Pharma, Berlin, Germany) using an automatic injector. Post-contrast images were obtained at 90, 180, 270, and 360 s after injection.
The whole tumor was semi-automatically segmented on the peak enhanced phase of CE-MRI by three-dimensional (3D) Slicer software (https://www.slicer.org/). The tissue surrounding or inside the tumor was also identified to provide information on both the microenvironment and the tumor itself. Specifically, the peritumoral area was obtained by expanding the tumor outward by a 2-pixel width (0.7 mm 9 2 pixels = 1.4 mm) and subtracting the tumor area, while the intratumoral area was obtained by shrinking by a 2-pixel width.
The extraction of spatial domain features, including shape features, first-order features, and texture features, was performed for each sample using the PyRadiomics package. Sequential features were calculated based on each spatial domain extracted from each CE phase (detailed descriptions of the radiomics features are presented in the electronic supplementary material). Feature selection was performed using the least absolute shrinkage and selection operator (Lasso) method (glmnet R package). Tuning parameter (k) selection in the Lasso model used tenfold cross-validation. The selected radiomics features were calculated for each patient using three different methods, including multivariate linear regression (glm R package) and support vector machine (SVM; e1071 R package). Radiomics prediction models were validated internally and externally. First, the trained classifiers were assessed by tenfold cross-validation via the glmnet R package. 26 The trained classifiers were then further tested in the validation cohort, and the area under the curve (AUC) of the receiver operating characteristics (ROC) curve was assessed.

Statistical Analysis
Student's t test and Wilcoxon's test were applied to compare continuous variables. Prior to the comparisons, the normality of variable distribution was tested using the Shapiro-Wilk test. Pearson's Chi-square test and Fisher's exact test were employed to compare unordered categorical variables. To explore the association between the HRD score and survival, Kaplan-Meier analysis and a Cox proportional hazards model were employed, and comparison of survival between groups was conducted via the logrank test. The predictive radiomics HRD model was validated using the SVM and logistic regression (LR) classifiers. First, the trained classifiers were assessed by tenfold cross-validation to reduce potential overfitting. The trained classifiers were then further tested in the validation cohort, and the AUC of the ROC curve was assessed. The ROC value was determined by the Youden Index. In multiple hypothesis testing, false discovery rate (FDR) correction was used to decrease false positive rates. All tests were two-sided, and a p-value \0.05 indicated significance unless otherwise stated. All of the statistical analyses were performed using R software (version 3.6.1; h ttp://www.R-project.org).

Study Population and Clinical Characteristics
A total of 343 TNBC patients with HRD scores were included in the present study, of whom 135 patients (39.4%) had a high HRD score and 208 (60.6%) had a low HRD score. The age and Ki-67 index parameters were found to differ significantly between the high and low HRD groups (p = 0.02 and 0.0004, respectively). No significant differences in pT or pN status were observed. The clinicopathological characteristics of the TNBC patients are listed in Table 1.

Prognostic Value of HRD in Triple-Negative Breast Cancer (TNBC)
Univariate Cox proportional hazards analysis showed that patients with high HRD scores had a higher probability of recurrence-free survival (RFS) in the overall TNBC cohort (p = 0.03) (Fig. 1a). Among the various TNBC subtypes, high HRD scores were correlated with a better RFS rate for the BLIS subtype only (p = 0.002) and not for the other subtypes (Fig. 1b). No significant correlations were observed between HRD and the prognosis of patients with the basal and nonbasal subtypes (Fig. 1c, d). A multivariate Cox proportional hazards model also revealed that HRD, together with pT and pN status, independently predicted RFS in all TNBC patients. Stage T3, stage N2/ N3, and low HRD scores were positively associated with worse RFS (Fig. 1e).
Correlation analysis indicated that cluster 3, characterizing immune-inflamed TNBC as revealed by our previous study, was more enriched in patients in the high HRD group (Fig. 2a). 27 Notably, among the distinct subtypes of TNBC, the LAR subtype was correlated with lower HRD scores, while the BLIS subtype was correlated with higher HRD scores; basal-like breast cancers were more enriched in the high HRD group (Fig. 2a). 7 In addition, patients in the high HRD group had a higher tumor mutation burden, more extensive lymphocyte infiltration, and fewer positive lymph nodes (Fig. 2b).

Transcriptomic and Metabolomic Characteristics of TNBCs with Low HRD Scores
Given the important role of the HRD score discussed in previous reports, we analyzed the transcriptomic differences between TBNC patients with low and high HRD scores (Fig. 3a). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that metabolic pathways were widely upregulated in patients with low HRD scores (Fig. 3b). Consistently, GSEA based on the KEGG and Reactome databases showed that metabolismrelated pathways, particularly amino acid and lipid metabolism pathways, were significantly upregulated, while DNA damage repair-and cell cycle-associated pathways were downregulated in patients with low HRD scores  (Fig. 3c, d). We then conducted a metabolomic analysis of samples with different HRD scores. Amino acids, peptides, and lipids were differentially expressed in patients with different HRD scores, which aligned with our findings from the transcriptomics analysis (Fig. 3e). Furthermore, lipidomics data revealed that sphingolipids were the largest proportion of lipids upregulated in patients with low HRD scores (Fig. 3f).

Predicting HRD Based on Noninvasive Radiomics
A total of 118 patients with preoperative MRI data included in the radiomic cohort were randomly assigned into training (n = 80) and validation cohorts (n = 38). Thirty-eight (32.2%) of these patients had high HRD scores and 80 (67.8%) had low HRD scores.
To investigate whether a noninvasive radiomics approach could be used to determine the HRD status with high accuracy, radiomic prediction models were curated via two steps: feature selection and model building. The HRD score, calculated as the arithmetic sum of three scores, was the gold standard to evaluate the prediction efficacy. 14 Three sequential radiomics features were selected to predict HRD status in the TNBC cohort (Fig. 4a), the mean value of the first-order minimum features after wavelet-LLL filtering from the tumor area, the kurtosis of the texture GLSZM Large Area High Gray Level Emphasis feature after wavelet-LHH filtering from the peritumoral region, and the difference in the first postenhanced phase and the plain scan value of the wavelet value of the first-order maximum features after wavelet-HHH filtering from the peritumoral region. The AUC, sensitivity, and specificity of SVM were 0.863, 0.797, and 0.857 in the training cohort and 0.739, 0.571, and 0.824 in the validation cohort, respectively, while the AUC, sensitivity, and specificity of LR were 0.864, 0.831, and 0.810 in the training cohort and 0.695, 0.571, and 0.882 in the validation cohort, respectively (Fig. 4b).
Age and Ki-67 were identified as variables for inclusion in the clinical model for predicting the HRD status, and achieved AUCs of 0.64 and 0.62 in the training and validation cohorts, respectively. The radiomics model had a significantly higher AUC than the clinical model (p \ 0.05), while the combined model showed no statistical significance in either cohort.

DISCUSSION
The HRD score is a useful index to evaluate HR and infer sensitivity to platinum chemotherapy and PARP inhibitors in TNBC patients. 14,18 By integrating the largescale multiomics data of TNBC patients treated at our institution, we elucidated the possible mechanisms underlying the distinct prognoses of patients with different HRD scores, and investigated the convenient clinical utilization of the HRD score in a noninvasive radiomics manner. Based on our multiomics database, we calculated the HRD score of each patient by arithmetically adding the three scores reported previously. 14 Patients in the high HRD group had significantly better survival outcomes than those with low scores. Furthermore, the use of the HRD score as an independent prognostic factor for TNBC was demonstrated by multivariable Cox regression. Transcriptomics analysis revealed that DNA damage repair-and cell cycleassociated pathways were upregulated in patients with high HRD scores, while metabolic pathways, particularly amino acid and lipid metabolism, were upregulated in patients with low HRD scores. We therefore investigated the metabolites upregulated in patients with low HRD scores using matched metabolomics data. The results showed that amino acids, peptides and lipids, particularly sphingolipids, were highly expressed in patients with low HRD scores, which provided a clue for targeting metabolites in TNBC patients with low HRD scores despite that they could not benefit from platinum chemotherapy or PARP inhibitors. Given the invasiveness and time required to determine HRD scores in the clinic, we established a radiomics signature to invasively predict the HRD status of TNBC patients.
HR is one mechanism of DNA DSB repair. 28 Tumors with higher HRD scores were considered to be more deficient in the process of HR, resulting in cell apoptosis. In the present study, we demonstrated that a higher HRD    29 Radiomics is an emerging field used to dissect medical imaging information based on high-throughput radiomics features. [30][31][32] In this study, we adopted noninvasive radiomics methods to differentiate TNBC patients with high and low HRD scores, thereby mitigating the cost and invasiveness of HRD detection. Previous studies have revealed the usefulness of radiomics in multiple clinical aspects. Aerts et al. revealed the prognostic value and correlation with gene expression patterns of noninvasive radiomics features in lung cancer and head and neck cancer. 33 The ability of radiomics to predict distant metastasis-free survival was assessed in locally advanced rectal cancer. 34 A 24-feature radiomics signature was constructed to predict lymph node status in colorectal cancer. A nomogram containing a radiomics signature performed well in predicting lymph node metastasis, with a C-index [0.75. 35 Notably, CD8-positive T cells and patient response to immunotherapy were predicted by radiomics features. 36,37 Our results, together with those of previous radiomics studies, provide support that radiomics has the potential to benefit clinical index prediction and decision making. In a previous study, our research team focused on the multiomics characteristics of TNBC patients in an Asian population. Transcriptomics analysis divided TNBC patients into four subtypes: LAR, BLIS, IM, and MES. 7 In the BLIS subtype, the HRD mutation signature was enriched, suggesting that platinum chemotherapy or PARP inhibitors for patients with high HRD scores, and more intensive therapies for patients with low HRD scores, are appropriate strategies for treating BLIS TNBCs with precision medicine. 12,13,[38][39][40][41][42] Based on the properties of each subtype, an umbrella trial was conducted, and 29% of refractory metastatic TNBC patients exhibited an objective response. 43 However, a consensus regarding treatment agents for TNBC patients with low HRD scores is lacking. Furthermore, our transcriptomics and metabolomics analyses demonstrated metabolic heterogeneity among TNBC patients. 44,45 These studies indicated the significance of transcriptome and metabolome information in the molecular landscapes and heterogeneity of TNBCs. Consequently, we conducted multiomics analyses of patients with different HRD statuses. We verified that cell cycle and DNA damage repair pathways were enriched in patients with high HRD scores and that amino acid and lipid metabolism pathways, particularly sphingolipid metabolism, were activated in patients with low HRD scores. These metabolic alterations might serve as druggable targets in TNBC patients with low HRD scores.
This study does have several limitations. First, our analysis was based on a retrospective cohort of a large multiomics dataset comprising transcriptomics, metabolomics, and radiomics data available; however, the results need to be retrospectively and prospectively validated in external cohorts and clinical trials. Second, the possible mechanisms underlying the worse prognosis of patients with low HRD scores were determined based on a multiomics cohort with clinicopathological information. Experimental investigations are necessary to further validate these findings and determine the explicit mechanisms.

CONCLUSIONS
We identified a set of patients with high HRD scores in the FUSCC multiomics TNBC cohort and demonstrated the prognostic value of the HRD score in TNBC. Notably, we proposed the practicability of predicting the HRD status with a noninvasive radiomics model. Moreover, amino acid and lipid metabolism, particularly sphingolipid metabolism, was shown to be enriched in the tumors of patients with low HRD scores, providing potential therapeutic targets for these patients.