3.1 Identification of differentially expressed genes
A total of 5079 DEGs were identified in 136 TNBC and 91 normal breast tissue samples (log2FC > = 1, Padj < 0.05). Figure 1B shows a volcano plot of the 2023 up-regulated and 3056 down-regulated genes. Figure 1C shows a heatmap of the top 50 DEGs with extreme values. A total of 199 significant differential NRGs were selected to construct a risk model in Fig. 1D.
3.2 Construction of a Prognostic Necroptosis-related genes Model
The univariate Cox regression analysis was conducted on 199 DEGs (Fig. 2A). The LASSO Cox regression analysis was employed to minimize the number of genes included in the prognostic signature, as depicted in Fig. 2B and Fig. 2C. A multivariate COX regression analysis was conducted to determine the regression coefficient of each gene. The model was constructed as follows: risk score = (-0.716 × expression level of GGT7) + (-1.644 × expression level of CCL25) + (− 0.779 × expression level of TNFRSF11B). All 136 TNBC samples were classified as high-risk (n = 68) or low-risk (n = 68) based on the median risk score.
3.3 Validation of the Prognostic Model in the TCGA cohort
The Kaplan-Meier curves from the TCGA cohort demonstrate a notable prognosis improvement (P < 0.05) in patients with low-risk ratings (Fig. 3A). The area under the curve (AUC) was determined through receiver operating characteristic (ROC) curve analysis. The AUCs for the risk model to predict overall survival at 3, 5, and 8 years were 0.891, 0.833, and 0.845, respectively (Fig. 3B). As illustrated in Fig. 3C, patients in distinct risk groups were significantly separated in two directions, as determined by t-SNE analyses. Figure 3D displays the rank distribution of risk scores, survival status, and expression patterns of the three NRGs in TNBC patients.
3.4 Validation of the Prognostic Model
The robustness of our predictive model was evaluated through external validation using the METABRIC cohort. The METABRIC dataset stratified patients into two groups based on their median TCGA scores: high-risk (n = 150) and low-risk (n = 144). The TCGA cohort identified a significant association between low-risk scores and favorable outcomes (Fig. 4A). According to Fig. 4B, the AUC of the ROCs is 0.602, 0.618, and 0.589 throughout 3, 5, and 8 years, respectively. As shown in Fig. 4C, the existence of two distinct subgroups of patients was confirmed by t-SNE analyses. Based on these findings, the predictive model was able to predict survival in this study accurately. According to Fig. 4D, the three NRGs in TNBC patients were ranked according to their risk scores, survival status, and expression patterns.
3.5 Independent Prognostic Value and Assessment
The independent predictive capability of the prognostic model was evaluated using univariate and multivariate Cox regression analyses. The efficacy of the prognostic model was evaluated using both the TCGA and METABRIC cohorts, which revealed that the risk score is an independent prognostic predictor in both cohorts (Fig. 5A and Fig. 5B). Furthermore, a nomogram was developed by integrating the three-gene predictive model with the stage to forecast the 3-year, 5-year, and 8-year OS in TNBC patients. The nomograms and calibration curves demonstrated the predictive model's utility in predicting outcomes (as illustrated in Fig. 5C-5D).
3.6 GO and KEGG Pathway Enrichment Analysis
To separate high-risk and low-risk groups, we annotated DEGs between high-risk and low-risk groups using GO and KEGG pathway enrichment analysis. Based on the analysis results, Fig. 6A and 6B present the top ten significantly enriched GO and KEGG terms. The GO enrichment analysis revealed that the DEGs were significantly enriched functional categories, including "negative regulation of immune system process," "extracellular structure organization," and "extracellular matrix organization.". Similarly, the KEGG pathway analysis confirmed the significant involvement of DEGs in pathways such as "Endocytosis," "Salmonella infection," "Focal adhesion," and "Cellular senescence," among others.
3.7 Analysis of the tumor immune microenvironment
To investigate the differences in immune cell infiltration between the two cohorts, an analysis was carried out to examine the relationship between the tumor immune microenvironment and the risk groups. The immune scores and ESTIMATE scores exhibited a significant difference between the high-risk and low-risk groups (P < 0.001) (Fig. 6A-6C). The low-risk group demonstrated a greater immune and ESTIMATE score than the high-risk group. Additionally, the low-risk group exhibited higher immune cell infiltration than the high-risk group.
We employed the CIBERSORT algorithm to analyze the infiltration of 22 immune cell types (Fig. 7D). The findings revealed that macrophage cells constituted the largest population among the 22 immune cell types, followed closely by resting T cells CD4 memory. Moreover, the expression of Macrophages M0 and Macrophages M2 cells was found to be higher in the high-risk group. Conversely, the low-risk group exhibited high expression levels of B cells naïve, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, and M1 cells. The study also examined the distribution of risk score, stromal score, immune score, and immune cell composition. In the low-risk group, immune regulatory pathways had a significant upregulation, such as antigen processing and presentation, TCR signaling, and BCR signaling, as illustrated in Fig. 7F. Furthermore, a correlation analysis was conducted between the risk score and 22 immune cells. The results demonstrated significant positive correlations between the risk score and eosinophils, Macrophages M0, Macrophages M2, activated mast cells, and monocytes. Conversely, we observed a negative correlation between the risk score and the frequencies of M1 and T follicular helper cells, as illustrated in Supplemental Fig. 1.
3.8 Analysis of immunotherapy
Previous analysis indicated significant differences between the two clusters in immune and stromal scores. Immune infiltration was significantly higher in the low-risk group than in the high-risk group. As immune cells play crucial roles in antitumor immune responses via mediating the lytic response, we quantified the expression of granzyme A (GZMA) and perforin (PRF1) to assess the variances in cytolytic activity between the two groups. The low-risk groups showed significantly higher cytolytic activity than the high-risk groups (Fig. 8A). We also determined the T cell inflammation scores (Fig. 8B) and found that the score was higher in the low-risk group than in the high-risk group. The dedifferentiation potential of tumor cells is evaluated via mRNAsi assessment. We compared the differences in mRNAsi (Fig. 8C) between the two groups, but no statistically significant results were detected.
Gene mutation data were analyzed to provide molecular insight into the immunological properties of high and low-risk groups. High-risk participants showed a statistically significant increase in mutations compared to the low-risk participants (as depicted in Fig. 9A-9B). A diverse range of gene mutations was identified, with missense mutations being the most prevalent, followed by nonsense mutations and frameshift deletions. The SNP was confirmed as the most fundamental variant type, while C > T was confirmed as the most fundamental SNV class. Additionally, we comprehensively analyzed the top ten highest mutation rates in both groups, demonstrating that TP53 and TTN had mutation rates exceeding 20%. The high-risk group exhibited frequent mutations in MUC16, PIK3CA, SPTA1, SYNE1, MUC17, RELN, BIRC6, and APOB, while the low-risk group frequently displayed mutations in CCDC168, MUC16, FAT3, DMD, PKHD1L1, USH2A, DYNC2H1, and MUC4 (Fig. 9C-9D).
To assess the efficacy of immunotherapy, we analyzed the correlation between tumor mutational burden (TMB) and risk score. Based on the positive correlation between risk score and TMB, the samples were classified into HRisk-HTMB, HRisk-LTMB&LRisk-HTMB, and LHisk-LTMB. KM analysis (Fig. 9E) indicated statistical differences between the three new groups, and the proportion of each group is depicted in Fig. 9F.
This study analyzed high and low-risk groups regarding the expression of seven common immune checkpoint genes. The results, as depicted in Boxplots (Fig. 10A-10B), revealed that individuals in the low-risk group exhibited significantly higher levels of TNFRSF9, CTLA4, PDCD1, TIGIT, IDO1, CD274, and LAG3 compared to those in the high-risk group. This finding suggests that patients in the low-risk group may not exhibit a favorable response to the targeted molecular immunotherapy. To further investigate the relationship between immune cell infiltrates and immunotherapy, the TIDE algorithm was employed to predict immunotherapy responses. The study findings indicate that a higher TIDE score may increase the probability of tumor immune evasion and decrease the effectiveness of ICI treatment. Our investigation revealed that the HRisk-HTMB group exhibited significantly lower TIDE scores compared to the LHisk-LTMB group (P < 0.05) (Fig. 10C), suggesting that the former group experienced greater benefits from ICI therapy. Furthermore, we observed a significantly higher T-cell dysfunction score (Fig. 10D) in the LHisk-LTMB group compared to the HRisk-HTMB group. Additionally, the Exclusion score (Fig. 10E) and MSI score (Fig. 10F) did not exhibit any significant differences
3.9 Prediction of drug sensitivity
After studying the content related to immunotherapy, we investigated the high-risk group of potential drugs with particular properties. After the screening, we selected the high and low-risk groups that exhibited elevated sensitivity to 4 inhibitor drugs for further investigation. First, we compared the IC50 values of four different drugs between the high-risk and low-risk groups. Subsequently, we listed and mapped the subtypes of all the specific drugs (Fig. 11). The IC50 scores showed a significant difference among the four inhibitors in both the high-risk and low-risk groups.
3.10 Validation of Gene Expression
We collected clinical specimens from TNBC patients and normal tissue to validate the expression profile of NRGs in TNBC patients. These specimens were analyzed using RT-PCR. Our results show that CCL25 was significantly upregulated in tumor samples, whereas there was no significant difference in expression levels of GGT7 and TNSRSF11B between tumor samples and control samples (Fig. 12A). In addition, the expression of each protein was verified by Western blot analysis (Fig. 12B). The Western blot assay revealed a higher increase in the expression of the CCL25 gene in tumor samples. On the other hand, there were no significant differences in the expression levels of GGT7 or TNSRSF11B between tumor specimens and controls (p > 0.05).