The signature of SARS-CoV-2-related genes predicts the immune therapeutic response and prognosis in breast cancer

Abstract

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an exceptionally contagious single-stranded RNA virus with a strong positive contagion. The COVID-19 pandemic refers to the swift worldwide dissemination of SARS-CoV-2 infection, which began in late 2019. The COVID-19 epidemic has disrupted many cancer treatments. A few reports indicate that the prevalence of SARS-CoV-2 has disrupted the treatment of breast cancer patients (BCs). However, the role of SARS-CoV-2 in the occurrence and prognosis of BC has not been elucidated. Here, we applied bioinformatics to construct a prognostic signature of SARS-CoV-2-related genes (SCRGs). Specifically, weighted gene co-expression network analysis (WGCNA) was utilized to extract co-expressed genes of differentially expressed genes (DEGs) in breast cancer and SCRGs. Then, we utilized the least absolute shrinkage and selection operator (LASSO) algorithm and univariate regression analysis to screen out three hub genes (DCTPP1, CLIP4 and ANO6) and constructed a risk score model. We further analyzed tumor immune invasion, HLA-related genes, immune checkpoint inhibitors (ICIs), and sensitivity to anticancer drugs in different SARS-CoV-2 related risk subgroups. In addition, we have developed a nomination map to expand clinical applicability. The results of our study indicate that BCs with a high-risk score are linked to negative outcomes, lower immune scores, and reduced responsiveness to anticancer medications. This suggests that the SARS-CoV-2 related signature could be used to guide prognosis assessment and treatment decisions for BCs.

Introduction

The COVID-19 pandemic has resulted in unparalleled societal turmoil and prompted swift changes in healthcare systems worldwide. Moris D, et al. propose that the COVID-19 pandemic has the potential to exacerbate pre-existing inequalities, namely for individuals with cancer [1]. All age groups are susceptible to infection. However, the severity of diseases is typically correlated with advanced age and pre-existing immune-compromising conditions, such as cancer [2]. Researches show that cancer patients have been particularly hard-hit by the impact of the COVID-19 pandemic. This impact encompasses adverse outcomes for cancer patients who contract COVID-19, disruptions to cancer treatment due to the pandemic, and significant setbacks for cancer research [3]. The dysregulated immune response in patients with COVID-19 has been a recurrent issue affecting symptoms and mortality rates. S100A8/A9 is primarily upregulated systemically by megakaryocytes and monocytes in peripheral blood, which may contribute to the cytokine storm frequently observed in severe patients [4]. Additionally, changes in immune characteristics targeting the respiratory tract are also key factors in the development of COVID-19. Studies have shown that pro-inflammatory monocyte-derived macrophages are abundantly present in the bronchoalveolar lavage fluid of severe COVID-19 patients, whereas moderate cases are characterized by the presence of highly clonally expanded CD8 + T cells [5]. Dysregulation of immune characteristics may also be another key factor in the progression of malignant tumors.

Breast cancer (BC), which continues to be the primary cause of cancer-related mortality among women and ranking among the top three most prevalent malignancies globally, is characterized by its high heterogeneity, resistance to treatment, and metastasis [6, 7]. For the majority of the last four decades, breast cancer incidence has increased [8].

The COVID-19 pandemic to medical services, including cancer care providers, has brought huge pressure [9]. Estimates based on good population surveillance models predict that the number of deaths related to breast cancer is likely to increase in proportion to the duration of cessation of screening [10, 11].

Although the aforementioned perspectives indicate significant disruptions to the treatment and diagnosis of BC due to the COVID-19 pandemic, there is a lack of reliable evidence regarding the involvement of SARS-CoV-2-related genes in the occurrence and prognosis of BC. Recent studies suggest that a high-risk score associated with SCRGs is linked to adverse outcomes in kidney renal clear cell carcinoma and glioma [12, 13]. Additionally, the research findings by Ren Q, et al. indicate that the risk signature associated with SCRGs plays a crucial role in predicting the prognosis of esophageal squamous cell cancer. A high-risk score signifies an unfavorable prognosis for patients [14]. However, the relationship between breast cancer and SARS-CoV-2 is yet to be elucidated. We speculate that SARS-CoV-2 may play an indispensable role in the occurrence of breast cancer and tumor-related immune regulation. In our study, we employed bioinformatics methods to construct a risk scoring model associated with SCRGs. We categorized 993 breast cancers (BCs) from The Cancer Genome Atlas (TCGA) database into two subgroups: high-risk and low-risk. Subsequently, we investigated differences between these two groups in terms of overall survival, immune cell infiltration, tumor immune dysfunction and exclusion (TIDE) score, as well as variations in the sensitivity to immune checkpoint inhibitors (ICIs) and anticancer drugs. Furthermore, we validated our risk scoring model in the GEO (GSE21653) dataset. To enhance the applicability for clinical use, we generated nomograms and corresponding calibration curves associated with prognostic risk factors. Taken together, our study suggests that the signature associated with SCRGs may serve as an effective tool for predicting the prognosis of breast cancer, responsiveness to immunotherapy, and sensitivity to anticancer drugs, which holds significant clinical relevance in guiding the treatment strategies for breast cancer.

Materials and methods

Clinical cohort and data preprocessing

Firstly, we incorporated clinical information from 993 BCs and 106 healthy individuals obtained from TCGA website (https://www.cancer.gov/ccg/research/genome-sequencing/tcga). Additionally, transcriptome and somatic mutation data were extracted for further analysis. Patients lacking clinical information were excluded from this study. The comprehensive summary of clinical baseline data is provided in Supplementary Table S1. Concerning the transcriptome data, the Level 3 HTSeq-FPKM format was converted to transcripts per million reads (TMP) format and log2(x + 1) transformed to mitigate expression differences among individuals. Furthermore, we included 266 BCs from the GEO database(https://www.ncbi.nlm.nih.gov/geo/) and subjected them to the same preprocessing methods as a validation set for further analysis. A collection of 333 genes associated with SARS-CoV-2 infection was compiled from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) and summarized in Supplementary Table S2. As the clinical information and transcriptome data designed for this study originate from publicly available databases, patient informed consent is not required.

Identification of DEGs related to SCRG in BCs

The “limma” package was utilized to identify differentially expressed genes (DEGs) between normal tissue and adjacent tissue in BCs (|log2FC| > 1, FDR < 0.05). Subsequently, a venn diagram was employed to capture 16 co-expressed genes shared between DEGs and SCRGs. Additionally, the GeneMANIA database (http://genemania.org/) was used to construct a protein-protein interaction network (PPI) for these 16 co-expressed genes. Interestingly, the “ggplot2” package was also applied to generate a correlation heatmap for the co-expressed genes.

WGCNA

Using the “WGCNA” package in R, we created a sample tree to identify outlier values and subsequently constructed a WGCNA of the BC expression matrix. The adjacency matrix was converted into a topology overlapping matrix (TOM). We used R2 = 0.90 as the criterion and applied Dynamic Tree Cut to generate modules, with parameters set to deepSplit = 2 and minModuleSize = 20. Finally, we calculated the correlation between different modules and the progression of breast cancer, identifying the modules with the strongest correlations (P < 0.05).

The construction of the risk signature and gene set variation analysis (GSVA)

The LASSO algorithm and univariate Cox regression analysis were employed to further screen prognosis-related genes from the 16 co-expressed genes. In addition, the “survival” package was also employed to screen prognostic markers among the 16 co-expressed genes. We intersected the genes selected through LASSO and univariate Cox regression analysis to determine the hub genes. It is noteworthy that we calculated the risk scores for each breast cancer patient by utilizing the coefficients obtained from the LASSO algorithm. The risk score for each individual is determined by multiplying the expression level of each hub gene by the sum of its LASSO coefficients: Riskscore = ∑Coefi x Exp_i. Subsequently, we divided the 993 breast cancer patients into high and low-risk subgroups using the median of the risk scores as the cutoff point. For the prognostic Kaplan-Meier curves, we employed the Log-rank algorithm to calculate the survival rates between the two distinct risk groups. Simultaneously, we utilized the “maxstat” package to compute the optimal cutoff value, ensuring the minimal p-value (with a grouping range set at 25-75%). Regarding GSVA, we utilized the “GSVA” and “limma” packages to compute differential signaling pathways between high and low-risk groups [15]. Subsequently, we employed the “ggplot2” package to create bar plots.

The correlation between the risk score and clinical features

Firstly, we used the “ggplot2” package to generate box plots for visualizing the differences in risk scores among different clinical subgroups of BCs. Next, we compared the overall survival (OS) between high and low-risk groups within different clinical subgroups. Furthermore, we utilized heatmaps in both the TCGA training set and GSE21653 validation set to visually represent the relationship between risk scores and prognosis.

Tumor immune analysis related to risk score

For tumor immune analysis between different risk subgroups, we employed the CIBERSORT algorithm and ESTIMATE algorithm to assess immune cell infiltration. Additionally, we analyzed the differences in immunophenoscore (IPS) from individuals of different risk groups using data from The Cancer Immunome Database (TCIA, https://tcia.at/home). Tumor immune and T cell exhaustion scores were obtained from the TIDE database (http://tide.dfci.harvard.edu/).

Analysis of sensitivity to anticancer drugs and somatic mutations

The “pRRophetic” package was applied to estimate the sensitivity of each breast cancer patient to chemotherapy drugs from the genomic of drug sensitivity in cancer database (GDSC; https://www.cancerrxgene.org/) [16]. Specifically, the half-maximal inhibitory concentration (IC50) was quantified using the “pRRophetic” package in R, utilizing IC50 to measure sensitivity to drugs. The differences in IC50 values for chemotherapy drugs between different risk groups were compared. Additionally, the “maftools” package was used to analyze somatic mutation data from different risk groups.

Construction of a prognostic-related risk factor model

The “survival” package was used for univariate and multivariate Cox regression analysis of different prognostic risk factors. Additionally, the “rms” package was employed to generate a nomogram for prognostic-related risk factors. Prognostic calibration curves for 1-, 3-, and 5- years were generated to assess the reliability of the model.

Statistics

All data analyses were conducted using R (v 4.3.1). The Mann-Whitney U test was employed for non-parametric numerical data, the Chi-square test for categorical variables, and logistic regression for continuous variables. The Log-rank test was used to compare survival curves between two groups. The Wilcoxon rank-sum test was utilized to compare differences between variables within groups. A significance level of p < 0.05 was considered indicative of a statistically significant difference.

Results

Extracting DEGs-related SARS-CoV-2 co-expressed genes and PPI construction

To provide a comprehensive overview of our research, we have created a flowchart illustrating the step-by-step process of our study (Fig. 1). To understand the SCRGs associated with the occurrence of invasive breast cancer, we first identified DEGs between normal tissue and tumor tissue using the “limma” package. We found 586 upregulated DEGs and 1104 downregulated DEGs. A volcano plot was generated to visualize the 1690 DEGs (Fig. 2A). Subsequently, we performed WGCNA on these DEGs and identified that 127 genes in the blue module were most strongly correlated with BC (Fig. 2B). Additionally, we analyzed the correlation between the blue module and BC. The scatter plot indicated a strong correlation (r = 0.58, p < 0.001). To visually observe the genes involved in both breast cancer and SARS-CoV2, we extracted 1690 DEGs, 127 BC-related module genes, and 333 SCRGs from the Genecard website. The Venn diagram showed that a total of 16 co-expression genes were selected (Fig. 2D). To further understand the interrelationships among these 16 co-expression genes, we used the GeneMANIA online website(http://genemania.org/) to generate the PPI diagram of the 16 co-expression genes (Fig. 2E). Moreover, a correlation heatmap was also used to visualize the relationships of the 16 co-expression genes (Fig. 2F). The results indicated a strong correlation among the 16 co-expression genes.

Fig. 1

Flowchart of the study

Fig. 2

Extraction of DEGs related to SARS-CoV-2 co-expression genes and construction of the PPI network. A Volcano plot of DEGs between breast cancer tumor tissues and adjacent normal tissues. B Heatmap of WGCNA module correlations related to DEGs. C Scatter plot of the module (blue) with the strongest correlation to breast cancer. D Venn diagram of DEGs, WGCNA, and SARS-CoV-2 genes. E GeneMANIA report for the 16 co-expression genes. F Heatmap of correlations among the 16 co-expression genes. *p < 0.05, **p < 0.01, ***p < 0.001

Selection of hub genes

To further investigate the core genes affecting BCs prognosis, we used LASSO analysis and univariate Cox regression analysis to screen the 16 co-expression genes. LASSO analysis indicated that 13 genes were associated with prognosis (CEP112, AASS, ANO6, CLIP4, FBLN5, ATP6AP1, DCTPP1, CHPF, MYCBP2, NAT14, SPART, SLC9A3R1, PRKAR2B) (Fig. 3A). According to the Partial Likelihood Deviance plot, the lambda value was set (Fig. 3B). Furthermore, we confirmed the prognostic genes using univariate analysis. The results indicated that DCTPP1(p < 0.001), ANO6(p = 0.023), and CLIP4(p = 0.037) were significantly associated with BCs prognosis (Fig. 3C). Univariate Cox regression analysis indicated that DCTPP1 and ANO6 were associated with poor prognosis of BCs, while CLIP4 was associated with favorable outcome of BCs. Subsequently, we selected the intersecting genes from LASSO and univariate Cox regression analyses (DCTPP1, ANO6, and CLIP4) (Fig. 3D), and calculated the risk score using the following formula: Riskscore = 0.347237913 * DCTPP1 + 0.406875531 * ANO6–0.144021618 * CLIP4. Finally, we divided all BCs into high and low-risk groups based on the median risk score. Furthermore, to explore the expression differences in signaling pathways between different risk score groups, we employed the GSVA algorithm to assess signaling pathways in all individuals. The results revealed that, compared to the low-risk group, individuals in the high-risk group exhibited upregulation in adipogenesis, PI3K/AKT/mTOR signaling, p53 pathway, oxidative phosphorylation, and Notch signaling, while hypoxia, inflammatory response, interferon alpha/gamma response were downregulated (Fig. 3E).

Fig. 3

Hub gene selection. A LASSO-Cox analysis coefficient plot. B Partial Likelihood Deviance plot of LASSO-Cox analysis. C Forest plot of univariate prognosis analysis for 16 co-expressed genes. D Venn diagram comparing genes from LASSO analysis and univariate analysis. E GSVA analysis bar plot depicting pathway correlations in high and low-risk groups

Correlation analysis between risk score and clinical features

To explore the relationship between risk scores and clinical features, we analyzed the association between risk scores and prognosis in different clinical subgroups. Firstly, our results indicated that risk scores were unrelated to tumor volume and distant metastasis status (Fig. 4A, C). However, BCs with lymph node metastasis and higher pathological grade showed higher risk scores compared to those without lymph node metastasis or with lower pathological grade (Fig. 4B, D). To further understand the impact of high and low-risk scores on BCs prognosis in different clinical subgroups, we generated prognostic Kaplan-Meier survival curves using the Log-rank algorithm. Our results suggested that in subgroups with different tumor volumes, no lymph node metastasis, no distant metastasis, and different pathological grades, the high-risk score group had a poorer prognosis compared to the low-risk score group (Fig. 4E-J). Additionally, we plotted prognostic survival curves for high and low-risk groups in the TCGA training set and GSE21653 validation set. The results were consistent with the previous findings, indicating an unfavorable prognosis for the high-risk score group (Fig. 4K, L).

Fig. 4

Correlation analysis between risk scores and clinical features. A-D Differential analysis of risk scores in subgroups based on different clinical features. Comparison of risk scores in subgroups based on tumor volume (A), lymph node metastasis (B), distant metastasis (C), and pathological grade (D). E-J Kaplan-Meier survival curves for high and low-risk subgroups under different clinical features. K-L Prognostic survival curves for high and low-risk groups of patients in TCGA training set (K) and GSE21653 validation set (L)

Prognostic riskscore and somatic mutation analysis

To further validate the relationship between risk scores and BCs prognosis, we generated prognostic risk score heatmaps from the TCGA training set and validated them in the GSE21653 dataset (Fig. 5A, B). Our results indicated that a high-risk score was associated with shorter survival time and fewer surviving patients. To investigate the relationship between risk scores and the somatic mutation burden of BCs, we analyzed somatic mutation data from different risk subgroups. The waterfall plot results suggested that the number of somatic mutations in the low-risk subgroup was higher than that in the high-risk subgroup. Regarding mutated genes, TP53, PIK3CA, TTN, CDH1, and GATA3 were the top 5 genes with the highest mutation frequencies, accounting for 39.3%, 35.8%, 22.0%, 16.0%, and 14.3% of the total mutations, respectively. Furthermore, GATA3 had the highest insertion mutation frequency (Fig. 5C).

Fig. 5

Prognostic risk factors and somatic mutation analysis. A Prognostic risk score heatmap from the TCGA training set. B Prognostic risk score heatmap from the GSE21653 validation set. C Waterfall plot of somatic mutations in different risk score groups

Tumor-associated immune analysis

To comprehensively understand the relationship between risk scores and the tumor immune microenvironment, we employed the CIBERSORT algorithm to assess 22 tumor-associated immune cells (Fig. 6A). The results indicated that the high-risk score group had higher levels of M2 macrophages, neutrophils, and resting mast cells compared to the low-risk group. Conversely, B cells (memory), plasma cells, CD8 + T cells, and T cells (follicular helper) showed higher CIBERSORT scores in the low-risk score group. To visualize the immunological landscape of these 22 immune cells between the two subgroups, we generated a correlation heatmap and stacked plot (Fig. 6B, C). Additionally, we applied the ESTIMATE algorithm to assess TCGA patients. The results showed that the immunescore and ESTIMATEscore were lower in the high-risk subgroup compared to the low-risk subgroup, while the stromalscore showed no significant difference between the two groups (Fig. 6D). These findings suggest that BCs in the high-risk group exhibit a lower immunological score, with a marked enrichment of immune cells promoting tumorigenesis and a notable deficiency of tumor-associated immune cells. This may be associated with the unfavorable prognosis observed in the high-risk group.

Fig. 6

Analysis of the relationship between risk scores and tumor-associated immune cells. A Scores of 22 tumor-associated immune cells calculated using the CIBERSORT algorithm in different risk score groups. B-C Correlation heatmap (B) and stacked plot (C) of the 22 immune cells. D Stromal score, immune score, and ESTIMATE score in different risk score subgroups. *p < 0.05, **p < 0.01, ***p < 0.001

Evaluation of immunotherapy responsiveness in different riskscre subgroups

To investigate the sensitivity of different risk subgroups to immunotherapy, we evaluated immunoinhibitors, immunomodulatory molecules MHC, immunophenoscore, and immunostimulators in BCs. The results indicated that the expression of most immunoinhibitors, such as BTLA, LAG3, and PDCD1, was upregulated in the high-risk group, while TGFBR1 expression was downregulated (Fig. 7A). Furthermore, the expression of HLA-related genes was lower in the high-risk group compared to the low-risk group (Fig. 7B). Regarding immunophenoscore (IPS), the high-risk group showed higher scores for IPS, PD1/PDL1 blocker, CTLA4 blocker, and PD1/PDL1 & CTLA4 blocker compared to the low-risk group (Fig. 7C). As for immunostimulators, our results suggested that almost all relevant genes were downregulated in the high-risk group, except for CD276 and TNFSF4 (Fig. 7D). Additionally, in TIDE analysis, we found that BCs with no response to immunotherapy had higher risk scores (Fig. 7E). This implies that the high-risk group exhibits a poorer response to immunotherapy. In terms of TIDE scores, the high-risk group is associated with higher TIDE, Exclusion, MSI, TAM M2 scores, and lower IFNG and Dysfunction scores (Fig. 7F).

Fig. 7

Assessment of immune therapy responsiveness in different risk subgroups of breast cancer (BCs). A-D Differential expression of various immunoinhibitors (A) immunomodulatory molecules MHC (B), Immunophenoscore (C), immunostimulators (D). E Responsiveness to immune therapy in different risk subgroups. F TIDE analysis results for different risk subgroups

Sensitivity analysis of different risk score subgroups to anticancer drugs

To explore the sensitivity of BCs in high and low-risk groups to anticancer drugs for clinical guidance, we used the “pRRophetic” package to assess drug sensitivity from the GDSC database. The results suggested that BCs in the high-risk group exhibited lower sensitivity to gemcitabine (Fig. 8A), doxorubicin (Fig. 8B), docetaxel (Fig. 8C), cisplatin (Fig. 8D), vinorelbine (Fig. 8E), paclitaxel (Fig. 8F), methotrexate (Fig. 8H), and camptothecin (Fig. 8I), while they showed higher sensitivity to MG132(Fig. 8G). The above results may be associated with the adverse prognosis in the high-risk group.

Fig. 8

Sensitivity analysis of anticancer drugs in different risk score subgroups. Sensitivity analysis of different risk score groups to Gemcitabine (A), Doxorubicin (B), Docetaxel (C), Cisplatin (D), Vinorelbine (E), Paclitaxel (F), MG132 (G), Methotrexate (H), Camptothecin (I). *p < 0.05, **p < 0.01, ***p < 0.001

Construction of the prognostic factor prediction model related to riskscore

To identify prognostic factors related to BCs, we conducted univariate and multivariate Cox regression analysis on clinical features. The univariate Cox regression results indicated that age, lymph node metastasis, distant metastasis, lymph node examined count, the number of lymph nodes positive by HE, and the risk score were associated with prognosis (Fig. 9A). Subsequently, the factors selected from the single-factor analysis were further analyzed using multivariate Cox regression, revealing that age, lymph node metastasis, distant metastasis, and the risk score were independent prognostic factors for BCs (Fig. 9B). After selecting the prognostic risk factors, we constructed a prognostic prediction model to assess BCs prognosis (Fig. 9C). Additionally, we generated nomograms to predict the 1, 3, and 5-year survival rates for BCs (Fig. 9D). The model calculates a total score based on the clinical features of each individual and predicts the patient’s survival rate. Furthermore, to assess the reliability of the model, we also constructed calibration curves for the 1, 3, and 5-year prognostic models. The closer the curve is to the ideal line, the higher the model’s fit and accuracy. The results suggest that the prognostic prediction model is reliable.

Fig. 9

Construction of the prognostic factor prediction model related to risk scores. A Univariate Cox regression analysis of prognostic risk factors. B Multivariate Cox regression analysis of prognostic risk factors. C Column chart of prognostic risk factors. D Calibration curve of the prognostic model for 1, 3, and 5 years

Validation of hub gene diagnosis and prognostic value

To comprehensively assess the role of the three hub genes in tumorigenesis, we evaluated the expression of DCTPP1 in pan-cancer transcriptome data obtained from the UCSC database(http://genome.ucsc.edu/). From the results, it can be observed that DCTPP1 is upregulated in most cancers, which suggests that it may be an independent risk factor for malignancies (Fig. 10A). In addition, the KM survival curve indicates that upregulation of DCTPP1 is associated with unfavorable outcomes in breast cancer (Fig. 10B). To investigate the prognostic and diagnostic value of hub genes in BC, we separately analyzed the prognostic survival curves and diagnostic receiver operating characteristic (ROC) curves for DCTPP1 in BCs. DCTPP1 has a good diagnostic value in breast cancer patients, with area under curve (AUC) values of 0.902 (Fig. 10C). These results suggest that hub genes may play a crucial role in the development of breast cancer. To further confirm the expression of hub genes in BC patients, we extracted the immunohistochemical analysis results of the three hub genes from the HPA database. Compared to normal breast tissue, the protein expression of DCTPP1 was increased in tumor tissue (Fig. 10D). Immunohistochemical results are consistent with the previous findings regarding the impact of hub genes on the prognosis of BCs.

Fig. 10

Validation of hub gene diagnosis and prognostic value. A Expression of DCTPP1 in pan-cancers. B KM survival curve (OS) of DCTPP1 in TCGA training set. C diagnostic ROC curve of DCTPP1. D Immunohistochemical results of DCTPP1 in normal breast tissue and breast cancer tissue in HPA database

Discussion

The outbreak of SARS-CoV-2 is a significant global public health event affecting human health since the 2020s. It was first reported in Wuhan, China, at the end of 2019 [17]. According to epidemiological assessments, approximately 80% of SARS-CoV-2 infections have a mild course, while the remaining 20% will progress to more severe conditions. The overall estimated mortality rate for COVID-19 is 2–3% [18]. Apart from increasing the risk of blood clotting, COVID-19 also imposes a psychological burden on patients [19, 20]. Additionally, immunocompromised populations, such as cancer patients or those with hematologic disorders, face an elevated risk of contracting and succumbing to COVID-19 [21]. With the continuous development and updating of COVID-19 vaccines, evidence suggests significant benefits for cancer patients, especially in protecting against severe COVID-19 [22].

To date, the mechanisms of viral immune evasion remain unclear, posing a significant challenge in developing effective clinical treatment strategies. As critical processes in transcriptional regulation, RNA processing mechanisms such as alternative polyadenylation (APA) and alternative splicing (AS) are crucial for regulating human genes in various infectious diseases. Some researchers suggest that SARS-CoV-2 proteins can bind to APA factors, affecting their gene expression levels and thus regulating APA. Moreover, APA may disrupt antigen presentation by the MHC in infected cells, further aiding SARS-CoV-2 in evading the host immune response [23]. Consequently, APA might be a better predictor than AS in COVID-19 patients. Additionally, reports indicate that the extent of splicing dysregulation in the host correlates with the severity of COVID-19 and affects drug-protein interactions. It is widely believed that virus infection-induced AS of host RNA can alter protein translation in host cells [24]. The CASA free platform developed by Chen et al. can display changes in host cell splicing patterns under microbial induction, helping users quickly identify splicing patterns and their potential regulators under specific conditions [25]. Besides APA and AS, N6-methyladenosine (m6A) also plays a crucial role in the functions and molecular mechanisms of SARS-CoV-2-infected cell lines. Studies have reported that m6A is associated with various clinical states in COVID-19 patients and may be involved in mechanisms of vaccine breakthrough infections [26].

Breast cancer is one of the most common malignant tumors worldwide, posing a major health concern. Drug resistance, recurrence, and metastasis are the three most critical challenges in the clinical prognosis of breast cancer patients [27,28,29,30]. It is well-known that breast cancer patients faced substantial challenges in treatment and management during the COVID-19 pandemic. However, our understanding of the impact of SARS-CoV-2 on breast cancer remains limited. Recently, there have been reports indicating that SCRGs play a crucial role in guiding the treatment and prognosis of various malignancies [31, 32]. Nevertheless, the relationship between SCRGs and the occurrence and prognosis of breast cancer has not been reported. In our study, we constructed, for the first time, a prognostic prediction model based on SCRGs’ model features. We analyzed the reactivity to immunotherapy in different risk groups, providing assistance for clinical treatment decisions for breast cancers. Specifically, using LASSO regression analysis and univariate regression analysis, we identified three prognostically relevant hub genes (ANO6, CLIP4, and DCTPP1) from 16 co-DEGs. We categorized all breast cancer patients into high- and low-risk groups based on calculated risk scores. The validation of the hub genes is consistent with previous reports. Specifically, Niu et al. demonstrated that knocking down miR-378a-3p levels can enhance the expression of DCTPP1, which can activate DNA repair signals in breast cancer to promote tumor proliferation [33]. Fan et al. suggested that methylated CLIP4 might be a novel prognostic and therapeutic biomarker for breast cancer [34]. Tang et al. indicated that ANO6 is downregulated in BC and is involved in influencing breast cancer prognosis through macrophage polarization [35]. Furthermore, we conducted clinical correlation analysis, somatic mutation analysis, gene set variation analysis, TIDE analysis, and anticancer drug sensitivity analysis for different subgroups. To enhance clinical applicability, we selected independent risk factors to construct prognostic prediction models for 1, 3, and 5 years.

The mechanism of breast cancer is complex, involving the activation and transduction of multiple signaling pathways, such as the PI3K/Akt/mTOR pathway [36, 37], Notch signaling [38, 39], p53 pathway [40]. GSVA results indicate that DEGs in the high-risk group of breast cancers are enriched in the PI3K/Akt/mTOR pathway, Notch signaling, and p53 pathway, implying an adverse outcome for the high-risk group. Additionally, in tumor immune analysis, TAM M2 cells, which promote tumor occurrence, play a crucial role in triple negative breast cancer (TNBC) patients [41, 42]. In this study, CIBERSORT analysis suggested that TAM M2 cells were higher in the high-risk group than in the low-risk group. According to ESTIMATE analysis, the immune score of the high-risk group was lower than that of the low-risk group. The ESTIMATE score refers to the Stromal score combined with the Immune score, which can represent the tumor purity of a certain malignant tumor [43]. Studies have shown that low tumor purity is associated with lower OS, which is consistent with the poor prognosis in the high-risk group in our study [44].

Immunotherapy is an important means of treating breast cancer, and TNBC is the subtype most suitable for immunotherapy [45]. With an increasing understanding of cancer cells evading the immune system, specific immune checkpoint inhibitors have emerged [46]. In recent years, drugs targeting the PD-1/PD-L1 axis have become a research focus in TNBC immunotherapy, with PD-1/PD-L1 inhibitors such as pembrolizumab and atezolizumab demonstrating good clinical efficacy. In the process of T cell activation, CTLA-4 is associated with the initial stage of T cell activity, while PD-1/PD-L1 is associated with the later stage of T cell activation [47, 48]. By comparing different subgroups of immune checkpoint inhibitor (ICI)-related genes, we found that the high-risk group had lower expression of PD-1/PD-L1 and CTLA-4 blockers than the low-risk group, indicating poor responsiveness to immunotherapy in the high-risk group. TIDE score is associated with tumor-related immune escape and tolerance. The higher the TIDE score, the greater the probability of individual immune escape [49]. It is noteworthy that we found higher TIDE scores in the high-risk group, and patients with higher risk scores had poorer responses to immunotherapy. These results suggest that high-risk group breast cancers may have poor responsiveness to immunotherapy. ICI and chemotherapy can synergize in multiple trials. Currently, some clinical studies are exploring the prospect of combining late-stage TNBC with immunotherapy and chemotherapy. Importantly, we found that the high-risk group had lower sensitivity to most chemotherapy drugs than the low-risk group, indicating a certain similarity in predicting the effectiveness of immunotherapy.

The mutation of the TP53 gene is closely associated with tumor immunity and serves as an effective biomarker for predicting the response of different types of tumors to immunotherapy [50]. Studies have indicated that compared to TP53 wild-type breast cancer patients, those with TP53 mutations exhibit significantly higher levels of immune infiltration, as well as increased activity in various immune cells, functions, and pathways. The heightened immune responsiveness may be linked to a better prognosis in TP53-mutated patients [51]. Interestingly, our results reveal that the frequency of TP53 somatic cell mutations is higher in the low-risk group than in the high-risk group, aligning with our aforementioned clinical outcomes.

In terms of clinical subgroup correlation, the high-risk group is associated with poorer clinical staging and grading. Additionally, through univariate Cox regression analysis and multivariate Cox regression analysis, we identified independent prognostic risk factors including age, N-stage, M-stage, and risk score. To enhance clinical applicability, we constructed prognostic nomograms and calibration curves to assess the prognosis of breast cancers with different clinical features. The calibration curves indicate that the model is reliable.

In summary, for the first time, we have constructed a model based on SCRGs features to assess the immunotherapy and prognosis of BCs. This model holds certain value for guiding clinical decisions in the treatment of BCs and evaluating treatment efficacy. Additionally, we identified three hub genes (ANO6, CLIP4, DCTPP1), providing a theoretical foundation for further exploration of how SARS-CoV-2 may influence the occurrence of BC. Despite the innovation in our study, there are some limitations. Firstly, our research primarily relies on bioinformatics analysis, lacking confirmation from in vivo and in vitro experiments, which is essential for a more comprehensive understanding of the molecular biology level and mechanisms involved. Furthermore, the analysis in this study involved all breast cancer patients without distinguishing specific subtypes, such as the comparison between TNBC and non-TNBC cohorts. Lastly, individual differences between BCs may impact the prognostic features associated with SCRGs, necessitating further validation in clinical cohorts to determine if these features can be applicable in clinical practice.

Conclusions

In our study, we have, for the first time, developed a prognostic risk prediction model related to SCRGs in BCs. Our results indicate that the high-risk score group of BCs is associated with adverse prognosis, lower immune scores, and low IPS scores. Additionally, concerning sensitivity to immunotherapy, the high-risk score group is associated with the downregulation of ICIs and HLA-related genes, lower tumor-related immune cell infiltration, and high TIDE score. Moreover, in terms of anticancer drugs, the high-risk group exhibits lower sensitivity compared to the low-risk group for the majority of anticancer medications. In summary, our research findings suggest that the SCRGs-related risk prediction model can guide clinical strategies and prognosis predictions for BCs.