Predicting the Clinical Outcome of Triple-Negative Breast Cancer Based on the Gene Expression Characteristics of Necroptosis and Different Molecular Subtypes

Abstract

Necroptosis, a kind of programmed necrotic cell apoptosis, is the gatekeeper for the host to defend against the invasion of pathogens. It helps to regulate different biological processes regarding human cancer. Nevertheless, studies that determine the impact of death on triple-negative breast cancer (TNBC) are scarce. Therefore, this paper has comprehensively examined the expression as well as clinical significance of necroptosis in TNBC. ConsensusClusterPlus was used to establish a stable molecular classification that used the expression regarding the necroptosis-linked genes. The clinical and immune characteristics of different subclasses were evaluated. Then, the weighted gene coexpression network analysis (WGCNA) assisted in determining key modules, and we selected the genes exhibiting obvious association with necroptosis prognosis through the relationship with prognosis. The univariate Cox regression analysis together with least absolute shrinkage and selection operator (LASSO) techniques served for the construction of the necroptosis-related prognostic risk score (NPRS) model, and the pathway characteristics of NPRS model grouping were further studied. Finally, the NPRS, taking into account the clinicopathological features, used the decision tree model for enhancing the prognostic model as well as the survival prediction. First, two stable molecular subtypes with different prognosis and immune characteristics were identified using necroptosis marker genes. Then, the key modules were identified, and 10 genes significantly related to the prognosis of necroptosis were selected. Then, the clinical prognostic model of NPRS was developed considering the prognosis-linked necroptosis genes. Finally, the NPRS model, taking into account the clinicopathological features, adopted the decision tree model for enhancing the prognostic model as well as the survival prediction. Herein, two new molecular subgroups considering necroptosis-linked genes are proposed, and an NPRS model composed of 10 genes is developed, which maybe assist in the personalized treatment and clinical treatment guidance of TNBC patients.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Breast cancer acts as a representative cancer type for females. Despite recent advancements in treatment, such as targeted therapy, hormone therapy, chemotherapy, and radiotherapy, breast cancer is still the primary cause of death related to cancer around the world [1]. High histological grade, proliferation rate, and ductal histology are characteristics of triple-negative breast cancer (TNBC). It is also linked to the deficient expression of the human epidermal growth factor receptor-2 (HER2), progesterone receptor, and estrogen receptor (ER) [2]. Approximately 10–20% of breast cancer patients have TNBC. Affected patients display high mortality and higher recurrence rates than those with any other breast cancer subtype [3]. The major systemic therapy for TNBC is routine chemotherapy. However, efforts to identify efficient targets to enhance clinical outcomes have been motivated by the dearth of molecular-focused based therapy and the dismal prognosis of TNBC patients.

Programmed necrosis, or necroptosis, is characterized by the permeability of the plasma membrane and modifications in the mitochondria, which cause the release of the cytoplasmic components into the extracellular environment and give rise to an inflammatory response [4]. Necroptosis refers to a death mechanism of necrotic cell independent of caspase. Its primary mediators through gene regulation are mixed-lineage kinase domain-like protein (MLKL), receptor-interacting protein kinase 1 (RIPK1), and RIPK3 [5]. There is mounting evidence that necroptosis is essential for controlling carcinogenesis and cancer progression. Najafov et al. noted that necroptosis can accelerate the metastasis of the tumors and T cell death [6]. By releasing the chemokine (C-X-C motif) ligand (CXCL), necroptosis in pancreatic cancer can encourage tumor cell motility and invasion [7]. The drug response of prostate cancer patients is influenced by necroptosis [4]. Necroptosis has been connected to antitumor immunity in cancer immunotherapy. According to Wang et al., the lncRNA SNHG1/miR-21-5p/TLR4 regulatory axis and the prognostic characteristics associated with necrosis were identified in gastric adenocarcinoma. According to Shen et al., the necroptosis signal is triggered in breast tissue cancer and is highly correlated with tumor growth and malignant tumor markers [8]. However, more information is required on the probable molecular mechanism of the TNBC and the exact function of necroptosis in the disease prognosis.

In this study, the mutational and immunological properties of stable molecular subgroups were identified by consensus clustering using genes relevant to necroptosis. Then, the WGCNA technique assisted in confirming the major modules, and such association was considered to select genes with close relation to necroptosis prognosis. The genes involved in necroptosis were discovered using univariate Cox regression analysis and the LASSO techniques. Then, a clinical prognostic model and risk model were developed, which could be useful in determining the prognosis of the TNBC patients and their particular treatment course.

2. Materials and Methods

2.1. Data Collection

All data regarding RNAseq, overall survival (OS) of the patients, and their TCGA-BRCA characteristics came from The Cancer Genome Atlas (TCGA) GDC API. The GSE103091 chip data set of the TNBC came from NCBI’s Gene Expression Omnibus (GEO) database, which yielded 238 samples. Furthermore, the data set containing the TNBC sequencing data (METABRIC, Nature 2012 & Nat Commun 2016) was collected from the cBioPortal website (http://http://www.METABRIC.org/), hereinafter referred to as the cBioPortal data set. 74 necroptosis-linked genes used in an earlier study were also used in this report [9].

2.2. Data Preprocessing

Processing of TCGA-BRCA data set is as follows: (1) normal samples were discarded; (2) samples of TNBC were retained according to the clinical information table; (3) the samples without the birth time and OS status were eliminated, and the samples with a survival $\text{time}<10$ years were retained. Processing of GSE103091 data set is as follows: (1) normal samples were discarded; (2) samples without their birth time and OS status were eliminated, and samples with an OS value of <10 years were retained; (3) the probe of chip data was converted into gene symbols; (4) the expression data of multiple genes corresponding to one probe were removed; (5) if a gene had multiple expression values, its average value was taken as the expression value. The following is the cBioPortal data set processing: (1) normal samples were removed; (2) samples of TNBC according to the clinical information table were retained; (3) the samples without the birth time and OS values were eliminated, and the samples with an OS value of <10 years were retained.

The TCGA and cBioPortal data sets (hereinafter referred to as RNAseq data sets) were combined, and the batch effect between different data sets was removed through the ComBat function of sva package [10].

2.3. Molecular Typing of the Necroptosis-Linked Genes

The ConsensusClusterPlus package served for consensus clustering for the building of a consistent matrix, thereby clustering and typing the samples [11]. We obtained samples’ molecular subtypes from the expression data regarding genes related to necroptosis. The Pam algorithm and “maximum” served for measuring the distance, carrying out 500 bootstraps. 80% of patients were included in the training set in each bootstrap process. We defined the cluster number in the range of 2-10, meanwhile calculating the consistent matrix and the cumulative distribution function (CDF) for obtaining the molecular subtype, thereby assessing the optimal classification.

2.4. Construction of WGCNA

The WGCNA process was used for constructing a gene coexpression network [12]. First, for developing the gene expression similarity matrix, the following equation was used to estimate the absolute value regarding the Pearson correlation coefficient between the $i$ and $j$ genes: ${\text{S}}_{ij}=1+\text{cor}{x}_i+y_i/2$. It has been noted that the coexpression network was in agreement with the scale-free network, that is to say, $\log k$ of a node with the connectivity, $k$, exhibited a negative link to $\log Pk$ of the node probability, while the $R^2$ value was >0.85. Thereafter, we converted the expression matrix to the adjacency matrix, which was then converted to the topological matrix. An average-linkage hierarchical clustering technique served for the clustering of genes related to the TOM. Based on the hybrid dynamic pruning tree, we defined a minimal no. of genes in every gene network module as 200. The eigengenes of each gene module were calculated after the dynamic cutting method had defined the gene modules. The modules were then clustered and analyzed, and the modules that were adjacent to one another were combined to form new modules. The Pearson correlation coefficient $b$ of each pair of genes is represented by a soft threshold parameter, $\beta$ [13]. This process further makes the strong relationship stronger, while making the weak relationship weaker at the index level: $a_{ij}=1+\text{cor}{x}_i+y_i{/2}^{\beta }$. The feature vector genes, or ME, serve as every module’s representative genes and define the overall gene expression levels of the modules as follows: $\text{ME}=\text{princomp}{x}_{ij}^q$ ($i$ denotes the gene in the modulus $q$ whereas $j$ indicates the microarray samples in the modulus $q$). The expression pattern of the gene in all the data samples and the ME expression profile of a particular vector gene was utilized to identify the gene in a module using the Pearson correlation. This was described as the module membership (MM), where ${\text{MM}}_{\text{i}}^q=\text{cor}{x}_i,{\text{ME}}^q$, while ME denoted the expression profile of the gene $i$.

2.5. Prediction of Immunotherapy Reactivity

The tumor immune dysfunction and exclusion (TIDE) algorithm assisted in verifying the predictive influence exerted by the immune checkpoint inhibitor score (IMS) on clinical reactivity exhibited by immune checkpoint inhibitors (ICIs). TIDE algorithm is a calculation technique predicting the responsiveness of the immune checkpoint blockade (ICB), based on the gene expression profile [14]. It focused on assessing 3 cell types which restrict T cell infiltration (M2 subtype of cancer-associated fibroblast (CAF), tumor-associated macrophages (TAM), and myeloid-derived suppressor cells (MDSCs)) and 2 different mechanisms describing the tumor immune escape (the score on the dysfunction of the tumor infiltrating cytotoxic T lymphocytes (CTLs) and the score on the rejection of the immunosuppressive factors on CTLs).

2.6. Gene Set Enrichment Analysis (GSEA)

GSEA assisted in the pathway analysis for examining the pathways involving different biological activities in molecular subtypes. GSEA employed all the candidate gene sets in the Hallmark database [15].

2.7. Invasion Abundance Calculation of TME Cells

Microenvironment Cell Populations-Counter (MCP-Counter) served for analyzing the scores of the 10 immune cells [16], while the single sample gene set enrichment analysis (ssGSEA) analyzed the scores of the 28 immune cells [17]. Meanwhile, the Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression (ESTIMATE) process served for evaluating the general immune microenvironment infiltration score [18].

2.8. Constructing the NPRS Scoring System to Evaluate TNBC Samples

Identification of molecular subtype-related modules. Here, the entire expression profile of RNAseq was used for WGCNA analysis, and the most relevant modules of molecular subtypes were identified as the “key modules.” Then, the genes present in the key modules were extracted, wherein the genes showing a significant prognosis were chosen as those associated with the necroptosis phenotype. The total gene number was decreased using the LASSO regression [19] and stepAIC [20], while the prognosis-linked genes were acquired. For the construction of the NPRS scoring system, the following formula was employed for calculating the NPRS score of every patient: $\text{NPRS}=\Sigma \beta \text{i}\times \text{Expi}$, where “$\text{Expi}$” represents necroptosis prognosis-associated genes’ expressions and “$\beta$” represents the Cox regression coefficient regarding corresponding gene. A threshold value of “0” was taken into account for dividing patients into group with low NPRS risk and group with high NPRS risk. The Kaplan-Meier (KM) technique assisted in constructing the OS curve for the prognosis analysis, and a log-rank test helped to determine the significant differences between the groups.

2.9. Cancer Stem Cell

The expression data regarding pluripotent stem cell (PSC) samples (embryonic stem cell (ESC) as well as induced PSC (iPSC)) from Progenitor Cell Biology Consortium (PCBC) database were used for predicting as well as calculating the stem cell index by one-class logistic regression (OCLR) method. Firstly, only the sample data of ESC and iPSC are kept, which are collectively referred to as SC samples. The Ensembl IDs of SC samples are converted into gene symbol, and only the genes encoding proteins are kept. We obtain 78 SC samples at last, and the expression profiles regarding 8087 mRNA genes in each sample are kept. For the obtained expression profile, the average value assisted in centralizing each sample. Finally, the OCLR method in R package GelNet (V1.2.1) was used for calculating the weight vector of each gene on the processed data.

3. Results

3.1. Molecular Typing Based on the Necroptosis-Linked Genes

First, the expression exhibited by genes associated with necroptosis was retrieved from the RNAseq expression profile matrix. Eight genes were found to have significant OS values in the low and high gene expression groups ($p<0.01$). The RNAseq data were then clustered using consensus clustering based on the eight necroptosis-linked genes related to prognosis, and the appropriate number of clusters was established using CDF. When cluster was chosen as 2, based on the CDF delta area curve, the clustering outcome presented a strong stability (Figures 1(a) and 1(b)). In the end, we selected $k=2$ for generating two molecular subtypes. After examining the prognostic traits, the results showed that these two molecular subtypes had notable prognosis differences. Furthermore, the C1 group exhibited a poor prognosis, and C2 showed a better prognosis while the individuals with the C1 subtype had a significantly greater death rate relative to the C2 subtype (Figure 1(c)). Additionally, using the same technique to molecularly type the data from the GSE103091, it was clear that these two forms of molecular typing had significantly different prognoses, which was in line with the training set (Figure 1(d)). When the survival status of the various subtypes in the two data sets was further analyzed, it was discovered that C1 had a greater death rate than C2, which was aligned with the C1’s poor prognosis (Figures 1(e) and 1(f)).

[figure(s) omitted; refer to PDF]

3.2. Immune Characteristics between Subtypes

We compared the immunological features exhibited by the two subtypes for better revealing the immune microenvironment differences between both these subtypes. Significant differences were noted among 22 different types of immune cells by using the ssGSEA function of GSEA to analyze the scores of 28 types of immune cells. Compared to C1, C2 showed a higher immunological score (Figure 2(a)). The scores of the 10 immune cells were then analyzed using the MCP-Counter, and 8 of them showed different scores, while the C2 subtype also presented high scores (Figure 2(b)). The results of ESTIMATE analysis, which was used to assess the overall immune microenvironment infiltration score, were in agreement with those of MCP-Counter and ssGSEA (Figures 2(c)–2(e)).

[figure(s) omitted; refer to PDF]

3.3. Immunotherapy/Chemotherapy Differences between Subtypes

Furthermore, we analyzed the difference in immunotherapy between the two subtypes. Firstly, the section focused on comparing the subtypes in terms of different immune checkpoint expressions. Obviously, both the subtypes exhibited a differential expression of the 14 immune checkpoint genes (ICGs) (Figure 3(a)). As presented in Figure 3(b), the differences in subtypes of immunotherapy were analyzed. The TIDE scores of the RNAseq queues did not show obvious difference for both the molecular subtypes. C2 showed a higher dysfunction score than C1, while C1 showed a higher exclusion score than C2. In addition, the response degree of subtypes in the RNAseq cohort to traditional chemotherapy drugs (sunitinib, paclitaxel, crizotinib, S-trityl-L-cystine, and CMK) was also analyzed. The results indicated that C2 was more sensitive to these five drugs (Figure 3(c)).

[figure(s) omitted; refer to PDF]

3.4. WGCNA Analysis Identifying Molecular Subtype-Related Gene Modules

RNAseq data set served for identifying molecular subtype-associated gene modules using the “WGCNA” R software package. First, we clustered the samples and selected the coexpression module (Figure 4(a)). For ensuring the network to be scale-free, we selected a $\beta =3$ value (Figures 4(b) and 4(c)). In addition, the following settings were used to merge similar clusters into the new modules, i.e., $\text{deepSplit}=2$, $\text{height}=0.25$, and $\text{minModuleSize}=200$. This yielded 10 modules in total (Figure 4(d)). It is worth noting that the grey module represented the gene set incapable of being aggregated with the other modules. Figure 4(e) highlights the gene statistics regarding every module. Further, the relationship between every module and the molecular subtype was analyzed, and obviously, the yellow module showed an obviously positive relationship with the C1 subtype and a significantly negative relationship with the C2 subtype (Figure 4(f)). The module membership (MM) and the gene significance (GS) of the genes included in the yellow module were significantly and positively related ($r=0.43$, $p<1\text{e}-5$) (Figure 4(g)). Through the function enrichment of genes in the yellow module based on the “ClusterProfiler” R software package, the top 10 pathways with the most significant enrichment were shown. The yellow module presented obvious enrichment in pathways like the neurotrophin signaling pathway, Huntington’s disease, and sphingolipid signaling pathway (Supplementary Figure S1). Finally, the yellow module which was highly related to molecular typing was identified as the key gene module related to molecular typing.

[figure(s) omitted; refer to PDF]

3.5. Determination of Genes Related to Necroptosis Phenotype

The RNAseq data were then divided by 7 : 3 with 0.7 being the proportion of the training data set for the yellow module genes that showed a strong correlation to the molecular subtypes determined by WGCNA in the previous step. The training data set that included these genes underwent univariate Cox regression analysis, finding 35 genes with a greater influence on prognosis ($p<0.01$). It had 13 “protective” genes and 22 “risk” genes (Figure 5(a)). Additionally, the 35 genes in the RNAseq data set were compressed using LASSO regression for lowering the gene no. in the risk model. First, we focused on examining the change trajectory presented by each independent variable. The no. of independent variable coefficients being close to 0 elevated as the value of lambda gradually increased (Figure 5(b)). We adopted10-fold cross-validation for constructing the model, meanwhile examining the confidence interval (CI) for each lambda (Figure 5(c)). The model was optimized at $\text{lambda}=0.0789$. Hence, we selected 20 genes as the target genes for the subsequent analyses when $\text{lambda}=0.0789$. Stepwise multivariate regression analysis and the Akaike information criteria (AIC) served for optimizing the model combining the results of the LASSO analysis of the 20 genes. Finally, 10 genes were identified as necroptosis-linked genes from LASSO analysis that affect prognosis (N6AMT1, ZNF79, PRICKLE3, DOCK6, EEFSEC, RPAP1, GDPD3, ABHD8, RAB11B, and FBXO33) (Figure 5(d)).

[figure(s) omitted; refer to PDF]

3.6. Clinical Prognosis Model Establishment and Validation

The risk coefficients of 10 genes were obtained (Figure 6(a)). Then, according to the formula defined by the score of necroptosis of our samples, the necroptosis-linked prognostic risk score (NPRS) of every sample was determined. In the RNAseq training, the RNAseq validation, and the entire RNAseq data sets, the RS of each sample was calculated, respectively, and the best cutoff was categorized into group with high risk and group with low risk, followed by building their KM and ROC curves, respectively. The results of the RNAseq training, RNAseq validation, and the entire RNAseq data sets revealed the poor prognosis of high NPRS and good prognosis of low NPRS, and the ROC curve had a high AUC (Figures 6(b)–6(d)). This result had also been verified in TCGA data sets, GSE103091, and cBioPortal (Supplementary Figure S2 A-C).

[figure(s) omitted; refer to PDF]

3.7. Pathway Characteristics between NPRS Groups

The section focused on analyzing the correlation between the biological functions and NPRS of different samples in the RNAseq data set. Accordingly, these pathways were negatively correlated with NPRS of samples, and these pathways were mainly tumor-related (KEGG_WNT_SIGNALING_PATHWAY, KEGG_TGF_BETA_ SIGNALING_PATHWAY) (Figure 7(a)). Next, we selected some important biological pathways to check for activation in different NPRS groups. It was noted that in the RNAseq queue, compared with NPRS-low, 5 pathways were activated and 22 pathways were inhibited in NPRS-high, and 6 pathways were activated and 14 pathways were inhibited in the GSE queue (Figure 7(b)). On the whole, KRAS signaling up and TNF-α signaling via NFκB, etc., were the main inhibited pathway in the NPRS-high group. Besides, we calculated the mRNAsi in high group and low data set and found that high group had higher mRNAsi (Figure S3).

[figure(s) omitted; refer to PDF]

3.8. NPRS in Combination with Clinicopathological Features for Enhancing the Prognostic Model and the Survival Prediction

The age, N stage, T stage, M stage, and the NPRS of patients with TNBC in the RNAseq cohort were considered for constructing the decision tree. As found, only NPRS, stage, and age could be still observed in the decision tree, and four different risk subgroups (i.e., C1, C2, C3, and C4) were determined (Figure 8(a)). Significant differences were noted in the OS values of the four risk subgroups, and the prognosis of the C4 subgroup was the worst (Figure 8(b)). We analyzed the distribution of the four risk subgroups, finding the risk subgroups C2, C3, and C4 included the NPRS-high patients, while the patients in the C1 group were NPRS-low patients (Figure 8(c)). In addition, some differences were also noted in the OS status of the patients in various risk subgroups, and the distribution of death cases in the C4 subgroup was the highest (Figure 8(d)). As revealed by the univariate and multivariate Cox regression analyses of the NPRS and clinicopathological traits, the NPRS was a prognostic factor with the largest significance (Figures 8(e) and 8(f)). As revealed by the univariate Cox regression analysis results, NPRS could result in poor survival rate of the patients ($p=0.0014$, $\text{HR}=11$, 95% CI: 2.5-48) and remarkably predict the poor OS (Figure 8(e)). The multivariate Cox regression analysis adopted the statistically significant characteristic factors in the univariate Cox regression, finding NPRS as an independent risk factor ($p=0.0022$, $\text{HR}=10$, 95% CI: 2.2-44) (Figure 8(f)). For quantifying the risk evaluation and the survival probability regarding TNBC patients, we integrated NPRS with other clinicopathological characteristics to establish a nomogram, as shown in Figure 8(g). From the model results, NPRS had the highest effect on the OS prediction. Further, the calibration curve served for assessing the prediction accuracy of model. Obviously, the 3 calibration points in 2, 3, and 5 years presented similar prediction calibration curve to standard graph (Figure 8(g)), which suggested that the nomogram displayed an effective prediction performance. Nomogram and NPRS had higher AUC (Figure 8(h)). Additionally, the decision curve analysis (DCA) served for assessing the model reliability. The nomograms also showed a significantly higher advantage than the other curves (Figure 8(i)). Compared with other single clinicopathological features, nomograms showed the strongest ability to predict survival (Figure 8(j)).

[figure(s) omitted; refer to PDF]

4. Discussion

Recently, the necroptosis signaling pathway has been considered as a significant event in the regulation of tumorigenesis and progression, which remarkably impacts tumor development, tumor necrosis, tumor metastasis, and tumor immune response, suggesting the potential of targeting necroptosis as a new tumor treatment [21]. It is reported that the current correlation between necroptosis and antitumor immunity is summarized: there are two strategies to trigger antitumor immunity through necroptosis: (1) inoculation of necrotic tumor cell vaccine: damage-associated molecular patterns (DAMPs) were released by necrotizing apoptotic tumor cells to promote the maturation of bone marrow-derived dendritic cells (BMDCs), cross excitation of the effector T cells, and the resulting cytotoxic effects. A higher concentration of interferon-γ was noted in this process, which may represent a different anticancer method used by the CD8⁺ T cells. (2) Inoculation of necrotic fibroblast vaccine: necrotic cells secrete the NF-κB-derived signals, which leads to DC cell activation, increases the antigen load of normal tumor cells, and increases the CD8⁺ T cell-mediated tumor control [22].

TNBC is mostly seen in young women. Its clinical course is invasive, the probability of visceral metastasis and brain metastasis is high, and it shows a worse prognosis compared to other forms of BC. In addition, the efficacy of endocrine therapy for hormone receptors and targeted therapy for blocking HER2 is poor. Chemotherapy is only the main treatment, but it is very susceptible to drug resistance, which seriously affects the prognosis. Immunotherapy can prolong the survival of patients. Potential biomarkers of TNBC immunotherapy response include high tumor mutation burden (TMB), TILs, and immune infiltration transcriptional characteristics [23]. In this study, two molecular subtypes (i.e., C1 and C2) of TNBC were defined by the expression profile of necroptosis-related genes. Significant prognostic differences were noted between the 2 molecular subtypes, and C1 showed a poor prognosis. The interplay between the many cell types that make up the TME is connected to carcinogenesis, tumor progression, treatment resistance, and immune infiltration matrix [24]. We determined the immune cell infiltration level in the RNAseq cohort of patients, finding higher immunological score of the C2 subtype relative to the C1 subtype. Cancer immunotherapy depending on the ICIs has shown great clinical success in recent years. The link between various antitumor immunity and cell death mechanisms has been recently discovered. Induced pyroptosis, necroptosis, and ferroptosis with ICI exhibit a synergistic increased antitumor effect [22]. After investigating the response of the subtypes to immunotherapy, most ICGs presented differential expression in the 2 molecular subtypes. Considering immunotherapy scores, the C2 group had a higher dysfunction score than the exclusion score of the C1 group. In addition, the C2 subtype was more sensitive to sunitinib, paclitaxel, crizotinib, S-trityl-L-cystine, and CMK. These findings suggest that patients with the C1 subtype who have poor immune cell infiltration should consider alternative treatment strategies because they have a bad prognosis and meanwhile show a poor response to pharmacological therapy.

Furthermore, LASSO regression and AIC algorithm were used to obtain the prognostic risk model constructed by 10 genes (N6AMT1, ZNF79, PRICLE3, DOCK6, EEFSEC, RPAP1, GDPD3, ABHD8, RAB11B, and FBXO33). DNA N6-methyladenine (6 mA) is seen to be a novel form of DNA methylation that is present in different eukaryotic cells. Studies have found that overexpression of methylase N6AMT1 increases the viability of hepatocellular carcinoma (HCC) cells, inhibits apoptosis, and enhances migration and invasion [25]. A different study found the N6AMT1 showed an obviously lower expression in the TNBC tissues relative to the normal tissues [26]. It was found that DOCK6 presented overexpression in the gastric cancer tissues, and the positive expression is related to gastric cancer metastasis, indicating that the gastric cancer patients showed a poor prognosis. Through miR-148b-3p, Rac1 and Cdc42 can be activated, which directly affect the motility of gastric cancer cells, and can serve as a new treatment strategy specific to gastric cancer [27]. EEFSEC presents an obvious upregulation in prostate cancer cells, and the high expression leads to poor prognosis, affecting the proliferation, migration, invasion, and cell cycle of 22Rv1 cells [28]. Naka et al. found that the deficiency of GDPD3 reduced the levels of some lysophosphatidic acid (LPA) and lipid mediators in chronic myeloid leukemia cells. Its deficiency also activated AKT/mTORC1 signaling pathway and cell cycle and inhibited the interaction of Foxo3a/β-catenin in the stem cell nucleus of chronic myeloid leukemia (PMID: 32943626). RAB11B-AS1 is a natural lncRNA upregulated in human BC, which enhances the BC cell in terms of the invasion and migration and promotes tumor angiogenesis and distant metastasis of breast cancer [29]. There are no relevant reports of ZNF79, PRICKLE3, RPAP1, ABHD8, and FBXO33 in cancer, and more research is still needed.

Recent research has demonstrated that during tumor formation, tumor necroptosis might be induced by death factors to engage in the RIPK1 and RIPK3 pathways [30]. A new study noted that RIPK1 is necessary for inflammatory disorders instead of for the growth or metastasis of tumors [31]. As a result, the fundamental mechanism behind tumor necroptosis during carcinogenesis is yet unknown. According to research by Baik et al., Z-DNA binding protein 1 (ZBP1), rather than RIPK1, is the primary regulator of tumor necroptosis in tumor formation in breast cancer. ZBP1 deletion prevents tumor cells from necroptosis during tumor formation, preventing tumor metastasis in the MVT-1 BC model. For determining the correlation between NPRS and biological functions, different tumor-related pathways (KEGG_WNT_SIGNALING_PATHWAY, KEGG_TGF_BETA_SIGNALING_PATHWAY) were negatively correlated with NPRS of the samples through ssGSEA analysis. Next, it was seen that the different NPRS groups activated different pathways. The results showed that in the RNAseq queue, compared with NPRS-low, 5 pathways were activated and 22 pathways were inhibited in the NPRS-high queue, and 6 pathways were activated and 14 pathways were inhibited in the GSE queue. On the whole, KRAS signaling up, TNF-α signaling via NF-κB, etc. are the main inhibited pathways in the NPRS-high group. KRAS signaling is associated with cancer progression in many cancers. A few recent studies indicated that the KRAS mutations are found in about 30% of TNBC patients [32]. Through cell-mediated apoptosis and immune processes, tumor necrosis factor-α (TNF-α) was a complex connector between inflammation and cancer. It has been stated that TNF-α is involved in cell proliferation, tumor migration, tumor metastasis, matrix degradation, invasion, and angiogenesis [33]. Many researchers who investigated the molecular mechanism of necroptosis highlighted the role played by the TNF signaling pathway. Usually, TNF activates proinflammatory genes through NF-κB signaling to induce an inflammatory response and activate downstream NF-κB and MAPK pathways, resulting in increased proinflammatory gene expression [34]. Therefore, necroptosis exhibits an obvious association with tumorigenesis and development, and its internal complex regulatory mechanism remains to be further studied.

Future researches shall be conducted to combine the fundamental experiments and be deepened from the functional perspective. The study failed to consider other factors due to the lack of basic data on clinical follow-up, especially the diagnostic details.

Two TNBC necroptosis-related subtypes were found in this investigation, and subtype verification and molecular identification were carried out using separate data sets. However, this study has some limitations, though different molecular subtypes were investigated and confirmed in two different data sets. These findings are based on retrospective research, so additional functional studies are required to investigate the molecular roles played by the biomarkers in the 2 molecular subtypes. To further elucidate the mechanism of these targets, it is also required for the validation by using tissue samples and clinical patients. Overall, this research offers several prospective biological targets for developing fresh immunotherapies, which might eventually support the individualized treatment of TNBC patients.

5. Conclusions

The prognosis of TNBC was accurately predicted using a new prognostic risk model that included 10 genes linked to necroptosis. Further investigation shall be conducted for fully understanding the intricate molecular functions of these 10 genes. This study further highlighted the connection between necroptosis-linked genes and TNBC prognosis and immunotherapy. The findings of this study might enable clinical TNBC patients to receive precise and individualized care.

Acknowledgments

This study was funded by the Zhejiang Provincial Natural Science Foundation of China (LQ20H160038) and by the Shenzhen Special Fund for Innovation and Entrepreneurship of Overseas High-Level Talents (KQJSCX20180330124428928).

References

[1] H. Y. Lin, P. Y. Chu, "Advances in understanding mitochondrial microRNAs (mitomiRs) on the pathogenesis of triple-negative breast cancer (TNBC)," Oxidative Medicine and Cellular Longevity, vol. 2021,DOI: 10.1155/2021/5517777, 2021.

[2] D. D. Singh, D. K. Yadav, "TNBC: potential targeting of multiple receptors for a therapeutic breakthrough, nanomedicine, and immunotherapy," Biomedicine, vol. 9 no. 8,DOI: 10.3390/biomedicines9080876, 2021.

[3] C. Vetrei, M. Passariello, G. Froechlich, R. Rapuano Lembo, E. Sasso, N. Zambrano, C. de Lorenzo, "Novel combinations of human immunomodulatory mAbs lacking cardiotoxic effects for therapy of TNBC," Cancers, vol. 14 no. 1,DOI: 10.3390/cancers14010121, 2022.

[4] G. L. Beretta, N. Zaffaroni, "Necroptosis and prostate cancer: molecular mechanisms and therapeutic potential," Cell, vol. 11 no. 7,DOI: 10.3390/cells11071221, 2022.

[5] J. Sprooten, P. De Wijngaert, I. Vanmeerbeek, S. Martin, P. Vangheluwe, S. Schlenner, D. V. Krysko, J. B. Parys, G. Bultynck, P. Vandenabeele, A. D. Garg, "Necroptosis in immuno-oncology and cancer immunotherapy," Cell, vol. 9 no. 8,DOI: 10.3390/cells9081823, 2020.

[6] A. Najafov, H. Chen, J. Yuan, "Necroptosis and cancer," Trends in Cancer, vol. 3 no. 4, pp. 294-301, DOI: 10.1016/j.trecan.2017.03.002, 2017.

[7] Y. Ando, K. Ohuchida, Y. Otsubo, S. Kibe, S. Takesue, T. Abe, C. Iwamoto, K. Shindo, T. Moriyama, K. Nakata, Y. Miyasaka, T. Ohtsuka, Y. Oda, M. Nakamura, "Necroptosis in pancreatic cancer promotes cancer cell migration and invasion by release of CXCL5," PLoS One, vol. 15 no. 1, article e0228015,DOI: 10.1371/journal.pone.0228015, 2020.

[8] F. Shen, X. Pan, M. Li, Y. Chen, Y. Jiang, J. He, "Pharmacological inhibition of necroptosis promotes human breast cancer cell proliferation and metastasis," OncoTargets and Therapy, vol. Volume 13, pp. 3165-3176, DOI: 10.2147/OTT.S246899, 2020.

[9] S. Xin, J. Mao, C. Duan, J. Wang, Y. Lu, J. Yang, J. Hu, X. Liu, W. Guan, T. Wang, S. Wang, J. Liu, W. Song, X. Song, "Identification and quantification of necroptosis landscape on therapy and prognosis in kidney renal clear cell carcinoma," Frontiers in Genetics, vol. 13, article 832046,DOI: 10.3389/fgene.2022.832046, 2022.

[10] X. Xie, E. C. Wang, D. Xu, X. Shu, Y. F. Zhao, D. Guo, W. Fu, L. Wang, "Bioinformatics analysis reveals the potential diagnostic biomarkers for abdominal aortic aneurysm," Frontiers in Cardiovascular Medicine, vol. 8, article 656263,DOI: 10.3389/fcvm.2021.656263, 2021.

[11] M. D. Wilkerson, D. N. Hayes, "ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking," Bioinformatics (Oxford, England), vol. 26 no. 12, pp. 1572-1573, DOI: 10.1093/bioinformatics/btq170, 2010.

[12] P. Langfelder, S. Horvath, "WGCNA: an R package for weighted correlation network analysis," BMC Bioinformatics, vol. 9 no. 1,DOI: 10.1186/1471-2105-9-559, 2008.

[13] S. Horvath, B. Zhang, M. Carlson, K. V. Lu, S. Zhu, R. M. Felciano, M. F. Laurance, W. Zhao, S. Qi, Z. Chen, Y. Lee, A. C. Scheck, L. M. Liau, H. Wu, D. H. Geschwind, P. G. Febbo, H. I. Kornblum, T. F. Cloughesy, S. F. Nelson, P. S. Mischel, "Analysis of oncogenic signaling networks in glioblastoma identifies ASPM as a molecular target," Proceedings of the National Academy of Sciences of the United States of America, vol. 103 no. 46, pp. 17402-17407, DOI: 10.1073/pnas.0608396103, 2006.

[14] P. Jiang, S. Gu, D. Pan, J. Fu, A. Sahu, X. Hu, Z. Li, N. Traugh, X. Bu, B. Li, J. Liu, G. J. Freeman, M. A. Brown, K. W. Wucherpfennig, X. S. Liu, "Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response," Nature Medicine, vol. 24 no. 10, pp. 1550-1558, DOI: 10.1038/s41591-018-0136-1, 2018.

[15] A. Liberzon, C. Birger, H. Thorvaldsdóttir, M. Ghandi, J. P. Mesirov, P. Tamayo, "The molecular signatures database hallmark gene set collection," Cell Systems, vol. 1 no. 6, pp. 417-425, DOI: 10.1016/j.cels.2015.12.004, 2015.

[16] E. Becht, N. A. Giraldo, L. Lacroix, B. Buttard, N. Elarouci, F. Petitprez, J. Selves, P. Laurent-Puig, C. Sautès-Fridman, W. H. Fridman, A. de Reyniès, "Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression," Genome Biology, vol. 17 no. 1,DOI: 10.1186/s13059-016-1070-5, 2016.

[17] P. Charoentong, F. Finotello, M. Angelova, C. Mayer, M. Efremova, D. Rieder, H. Hackl, Z. Trajanoski, "Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade," Cell Reports, vol. 18 no. 1, pp. 248-262, DOI: 10.1016/j.celrep.2016.12.019, 2017.

[18] K. Yoshihara, M. Shahmoradgoli, E. Martínez, R. Vegesna, H. Kim, W. Torres-Garcia, V. Treviño, H. Shen, P. W. Laird, D. A. Levine, S. L. Carter, G. Getz, K. Stemke-Hale, G. B. Mills, R. G. W. Verhaak, "Inferring tumour purity and stromal and immune cell admixture from expression data," Nature Communications, vol. 4 no. 1,DOI: 10.1038/ncomms3612, 2013.

[19] J. Friedman, T. Hastie, R. Tibshirani, "Regularization paths for generalized linear models via coordinate descent," Journal of Statistical Software, vol. 33 no. 1,DOI: 10.18637/jss.v033.i01, 2010.

[20] Z. Zhang, "Variable selection with stepwise and best subset approaches," Annals of Translational Medicine, vol. 4 no. 7,DOI: 10.21037/atm.2016.03.35, 2016.

[21] Y. Gong, Z. Fan, G. Luo, C. Yang, Q. Huang, K. Fan, H. Cheng, K. Jin, Q. Ni, X. Yu, C. Liu, "The role of necroptosis in cancer biology and therapy," Molecular Cancer, vol. 18 no. 1,DOI: 10.1186/s12943-019-1029-8, 2019.

[22] R. Tang, J. Xu, B. Zhang, J. Liu, C. Liang, J. Hua, Q. Meng, X. Yu, S. Shi, "Ferroptosis, necroptosis, and pyroptosis in anticancer immunity," Journal of Hematology & Oncology, vol. 13 no. 1,DOI: 10.1186/s13045-020-00946-7, 2020.

[23] T. E. Keenan, S. M. Tolaney, "Role of immunotherapy in triple-negative breast cancer," Journal of the National Comprehensive Cancer Network: JNCCN, vol. 18 no. 4, pp. 479-489, DOI: 10.6004/jnccn.2020.7554, 2020.

[24] D. Hanahan, R. A. Weinberg, "Hallmarks of cancer: the next generation," Cell, vol. 144 no. 5, pp. 646-674, DOI: 10.1016/j.cell.2011.02.013, 2011.

[25] Q. Lin, J. W. Chen, H. Yin, M. A. Li, C. R. Zhou, T. F. Hao, T. Pan, C. Wu, Z. R. Li, D. Zhu, H. F. Wang, M. S. Huang, "DNA N6-methyladenine involvement and regulation of hepatocellular carcinoma development," Genomics, vol. 114 no. 2, article 110265,DOI: 10.1016/j.ygeno.2022.01.002, 2022.

[26] X. Sheng, J. Wang, Y. Guo, J. Zhang, J. Luo, "DNA N6-methyladenine (6mA) modification regulates drug resistance in triple negative breast cancer," Frontiers in Oncology, vol. 10, article 616098,DOI: 10.3389/fonc.2020.616098, 2021.

[27] X. Li, M. Jiang, D. Chen, B. Xu, R. Wang, Y. Chu, W. Wang, L. Zhou, Z. Lei, Y. Nie, D. Fan, Y. Shang, K. Wu, J. Liang, "miR-148b-3p inhibits gastric cancer metastasis by inhibiting the Dock6/Rac1/Cdc42 axis," Journal of Experimental & Clinical Cancer Research: CR, vol. 37 no. 1,DOI: 10.1186/s13046-018-0729-z, 2018.

[28] B. Xu, J. Hao, Q. Xie, N. Sa, S. Wang, H. Lu, P. Gao, G. Shi, X. Dong, "EEFSEC knockdown inhibits proliferation, migration and invasion of prostate cancer cells in vitro," Journal of Southern Medical University, vol. 41 no. 12, pp. 1787-1794, DOI: 10.12122/j.issn.1673-4254.2021.12.05, 2021.

[29] Y. Niu, L. Bao, Y. Chen, C. Wang, M. Luo, B. Zhang, M. Zhou, J. E. Wang, Y. V. Fang, A. Kumar, C. Xing, Y. Wang, W. Luo, "HIF2-induced long noncoding RNA RAB11B-AS1 promotes hypoxia-mediated angiogenesis and breast cancer metastasis," Cancer research, vol. 80 no. 5, pp. 964-975, DOI: 10.1158/0008-5472.CAN-19-1532, 2020.

[30] L. Seifert, G. Werba, S. Tiwari, N. N. Giao Ly, S. Alothman, D. Alqunaibit, A. Avanzi, R. Barilla, D. Daley, S. H. Greco, A. Torres-Hernandez, M. Pergamo, A. Ochi, C. P. Zambirinis, M. Pansari, M. Rendon, D. Tippens, M. Hundeyin, V. R. Mani, C. Hajdu, D. Engle, G. Miller, "The necrosome promotes pancreatic oncogenesis via CXCL1 and Mincle-induced immune suppression," Nature, vol. 532 no. 7598, pp. 245-249, DOI: 10.1038/nature17403, 2016.

[31] S. Patel, J. D. Webster, E. Varfolomeev, Y. C. Kwon, J. H. Cheng, J. Zhang, D. L. Dugger, K. E. Wickliffe, A. Maltzman, S. Sujatha-Bhaskar, P. Bir Kohli, S. Ramaswamy, G. Deshmukh, B. M. Liederer, R. Fong, G. Hamilton, P. Lupardus, P. Caplazi, W. P. Lee, M. van Lookeren Campagne, A. Johnson, B. S. McKenzie, M. R. Junttila, K. Newton, D. Vucic, "RIP1 inhibition blocks inflammatory diseases but not tumor growth or metastases," Cell Death and Differentiation, vol. 27 no. 1, pp. 161-175, DOI: 10.1038/s41418-019-0347-0, 2020.

[32] M. Kriegsmann, V. Endris, T. Wolf, N. Pfarr, A. Stenzinger, S. Loibl, C. Denkert, A. Schneeweiss, J. Budczies, P. Sinn, W. Weichert, "Mutational profiles in triple-negative breast cancer defined by ultradeep multigene sequencing show high rates of PI3K pathway alterations and clinically relevant entity subgroup specific differences," Oncotarget, vol. 5 no. 20, pp. 9952-9965, DOI: 10.18632/oncotarget.2481, 2014.

[33] D. Cruceriu, O. Baldasici, O. Balacescu, I. Berindan-Neagoe, "The dual role of tumor necrosis factor-alpha (TNF- α ) in breast cancer: molecular insights and therapeutic approaches," Cellular Oncology (Dordrecht), vol. 43 no. 1,DOI: 10.1007/s13402-019-00489-1, 2020.

[34] D. Bertheloot, E. Latz, B. S. Franklin, "Necroptosis, pyroptosis and apoptosis: an intricate game of cell death," Cellular & Molecular Immunology, vol. 18 no. 5, pp. 1106-1121, DOI: 10.1038/s41423-020-00630-3, 2021.