Necroptosis-related lncRNA to establish novel prognostic signature and predict the immunotherapy response in breast cancer.

BACKGROUND: Necroptosis is a type of programmed cell death, and recent researches have showed that lncRNAs could regulate the process of necroptosis in multiple cancers. We tried to screen necroptosis-related lncRNAs and investigate the immune landscape in breast cancer (BC).
METHODS: The samples of breast normal and cancer tissue were acquired from TCGA and GTEx databases. A risk prognostic model was constructed based on the identified necroptosis-related lncRNAs by Cox regression and least absolute shrinkage and selection operator (LASSO) method. Moreover, the forecast performance of this model was verified and accredited by synthetic approach. Subsequently, an accurate nomogram was constructed to predict the prognosis of BC patients. The biological differences were investigated through GO, GSEA, and immune analysis. The immunotherapy response was estimated through tumor mutation burden (TMB) and tumor immune dysfunction and exclusion (TIDE) score.
RESULTS: A total of 251 necroptosis-related lncRNAs were identified by differential coexpression analysis, and SH3BP5-AS1, AC012073.1, AC120114.1, LINC00377, AL133467.1, AC036108.3, and AC020663.2 were involved in the risk model, which had an excellent concordance with the prediction. The pathway analyses showed that immune-related pathways were relevant to the necroptosis-related lncRNAs risk model. And the risk score was significantly correlated with immune cell infiltration, as well as the ESTIMATE score. Most notably, the patients of higher risk score were characterized with increased TMB and decreased TIDE score, indicating that these patients showed better immune checkpoint blockade response.
CONCLUSION: These findings were conducive to understand the function of necroptosis-related lncRNAs in BC and provide a potential promising therapeutic strategy for BC.

INTRODUCTION

Worldwide, breast cancer (BC) is one of the most prevalent types of malignancy and a major cause of cancer death. The treatment landscape of patients with breast cancer has been rapidly evolving in recent years, and optimal therapy paradigm for breast cancer depends on subpopulations of patients.

Necroptosis, mainly mediated by receptor‐interacting protein kinase 1 (RIPK1), RIPK3, and mixed lineage kinase domain like pseudokinase (MLKL), belongs to the category of programmed cell death. , ,

Increasing evidence demonstrated that necroptosis played a pivotal role in the occurrence and progression of multifarious diseases, such as neurodegenerative diseases, ischemic cardiovascular, and cancer metastasis. Besides, necroptosis had dual impact on promoting and suppressing tumor growth in multiple tumor types. , , Accordingly, regulating tumor necroptosis may be a novel and potential therapeutic strategy for BC.

Long noncoding RNAs (lncRNAs) have attracted growing focus as tumor markers for early‐stage detection, diagnosis, prognosis, and prediction of drug therapy response. , Different lncRNAs regulate the expression of genes through epigenetic regulation or transcriptional alternation, and the aberrant expression of lncRNAs is closely linked to tumorigenesis, metastasis, and chemoresistance of cancer. , , Up to date, it has been reported that some lncRNAs engaged in modulating the process of necroptosis by interacting with miRNA to regulate necroptosis‐related genes products. , , Given that, further insight into the function of necroptosis‐related lncRNAs in BC may provide novel approach for precise treatment and individualized management.

MATERIALS AND METHODS

Normal and tumor sample extraction from dataset

The transcriptome RNA‐seq datasets (HTSeq‐Counts and FPKM) of female breast cancer (BC) and normal samples were acquired from The Cancer Genome Atlas (TCGA) and Genotype‐Tissue Expression Project (GTEx), respectively. The HTSeq‐Counts value matrix was used to search for the differentially expressed (DE) lncRNAs, while the FPKM values were transformed to TPM values for other analyses. After excluding the sample of male, 191 normal tissue samples (79 samples from GTEx dataset) and 1086 BC samples were obtained from two datasets. After ruling out the missing overall survival time, 1057 cases with survival time and 908 cases with full clinical pathology information were extracted for following analyses.

Identification of Necroptosis‐related lncRNA

According to necroptosis gene set M24779.gmt and previous literature search, 67 necroptosis‐related genes were collected for next identification. Then, we used GENCODE annotation file to identify 14,106 lncRNAs in the TCGA combined with GTEx datasets (http://cancergenome.nih.gov/abouttcga, http://www.gtexportal.org). The DE lncRNAs were screened by DESeq2 R package with the standard of |log2 fold change (FC) | >1, false discovery rate (FDR) <0.05, and p adjusted <0.05. The Pearson correlation algorithm was used to identify necroptosis‐related DE lncRNA with correlation filter >0.4 and p < 0.001. After the completion of screening steps, 251 necroptosis‐related lncRNAs were retrieved for further analysis. The analyses were based on limma R package.

Establishment and validation of prognostic risk assessment model

Preliminarily, the prognostic‐classified lncRNAs were selected by using univariate Cox (uni‐Cox) regression with p value < 0.05. Then, LASSO regression analysis was made to filter necroptosis‐related lncRNA with 10‐fold cross‐validation. Further, the necroptosis‐related lncRNAs screened by LASSO method were used for multivariate (multi‐Cox) proportional hazards regression and risk model construction. The risk score was calculated by using following formula: risk score =  expression (lncRNAk) × coefficient(lncRNAk). The median risk score that calculated by the above formula was used to stratify the BC patients into low‐ and high‐risk groups. The chi‐square test was used to validate the correlation of the clinical features and the risk group. The independent variables were assessed by uni‐Cox and multi‐Cox regression analyses, respectively. Receiver operating characteristic (ROC) curves and concordance index(C‐index) were subsequently applied to measure the precision of the model. The analyses were based on survival, caret, glmnet, rms, survminer, and timeROC R packages.

Predictive nomogram construction and calibration

With rms R package, the nomogram was set up based on risk group, age, and clinicopathological factors. The nomogram aimed to evaluating the predictive efficacy of risk score we got for 1‐, 3‐, 5‐, and 10‐year overall survival rates. Subsequently, the calibration curve was developed to illustrate the prediction power of the established nomogram model.

PCA, GO, and GSEA analysis

The principal component analysis (PCA) was used to classify BC samples through necroptosis‐related lncRNAs expression patterns. Additionally, the spatial distribution of samples was displayed through three‐dimensional scatterplot. We identified the differential genes among the low‐ and high‐risk groups for subsequent Gene Ontology (GO) analysis, aiming to investigate the relevant biological process. Furtherly, differentially expressed KEGG pathways in two groups were identified by Gene Set Enrichment Analysis (GSEA). The KEGG gene set(c2.cp.kegg.v7.0.symbols.gmt) was derived from the website (https://www.gsea‐msigdb.org/). The limma, org. Hs.eg.db, clusterProfiler, and enrichplot package based on R 4.1.1 were used for the analysis. The threshold of significantly enriched biological processes and pathways was set as p < 0.05 and FDR <0.25.

Tumor microenvironment in low‐ and high‐risk groups

The correlations between risk score and tumor‐infiltrating immune cells (TIICs) were evaluated by CIBERSORT algorithm. Furthermore, with ESTIMATE R package, we calculated the tumor microenvironment (TME) score, including stromal score, immune score, and estimate score between low‐ and high‐risk groups.

Tumor mutation burden and Tumor Immune Dysfunction and Exclusion score

The somatic mutation file (TCGA.BRCA.varscan.DR‐10.0.somatic) was obtained from the TCGA website. The original mutation annotation format (MAF) was divided into two groups according to the risk score. Then, we calculated the tumor mutation burden (TMB) score according to the somatic mutation data for each patient in the two groups. The foregoing analyses were based on maftools R package. Potential immune checkpoint blockade (ICB) response was assessed by tumor immune dysfunction and exclusion (TIDE) algorithm. Finally, we used pRRophetic R package to calculate the semi‐inhibitory concentration (IC50) values of chemotherapeutic drugs.

RESULTS

Altered Expression of the necroptosis‐related lncRNA in BRCA

We analyzed the necroptosis‐related lncRNA expression level in 1,086 BC samples and 191 normal samples from the TCGA and GTEx datasets, and 2,848 DE lncRNAs were obtained. Then, 251 necroptosis‐related lncRNAs were identified in the DE lncRNAs by the Pearson correlation algorithm.

Construction and verification of prognosis risk assessment model

Firstly, 16 lncRNAs was extracted by means of the uni‐Cox regression analysis. Then, 12 lncRNAs were acquired by LASSO analysis, and 7 of which were brought in the multi‐Cox proportional hazards model (Figure 1A‐E). Finally, we got the risk score with the formula from multivariate Cox regression: Risk score = SH3BP5‐AS1 × (−0.3537) + AC012073.1 × (0.3945) + AC120114.1 × (0.3010) + LINC00377 × (−1.6837) + AL133467.1 × (−0.7597) + AC036108.3 × (−0.3151) + AC020663.2 × (−0.5513). In the complete set, training and validation partitions, all patients in the high‐risk group had a significantly shorter overall survival duration (Figure 2A‐I). The same results were displayed in the different clinicopathologic characteristics (Figure 2M). The area under the 1‐,3‐,5‐, and 10‐year ROC curve (AUC) was 0.731, 0.643, 0.641, and 0.694, respectively (Figure 3A). At the 10‐year ROC of the model, the AUC of risk score was 0.731, demonstrating strong predictive ability compared with other clinicopathology features (Figure 3B). The 1‐year C‐index in the risk model was 0.726 (Figure 3C). In uni‐Cox and multi‐Cox regression, the hazard ratio (HR) of the risk score were 1.246 and 1.279, respectively (both p value < 0.001) (Figure 3D‐E).

FIGURE 1

Identification of prognostic necroptosis‐related lncRNAs in BC. (A) The prognostic lncRNAs extracted by uni‐Cox regression analysis. (B) The heatmap of 16 lncRNAs expression. (C) The 10‐fold cross‐validation for tuning parameter selection in the LASSO model. (D) The LASSO coefficient profile of 16 survival‐associated lncRNAs. (E) 7 lncRNAs identified by multi‐Cox proportional hazard model. (F) Correlations between lncRNAs in the risk model and necroptosis‐related genes

FIGURE 2

Prognosis of the risk model in the different sets. (A–C) Demonstration of risk model of the training, validation, and complete sets. (D–F) Survival time and clinical endpoint in the training, validation, and complete sets. (G–I) The heatmap of 7 lncRNAs expression in the training, validation, and complete sets. (J–L) K‐M survival curves of OS of patients between the two groups in the training, validation, and complete sets. (M) K‐M survival curves of OS prognostic value stratified by clinicopathologic characteristics in the complete set

FIGURE 3

Verification of prognosis risk assessment model. (A) The 1‐, 3‐, 5‐, and 10‐year ROC curves of the complete sets. (B) The ROC curves of risk score and clinicopathologic features. (C) The C‐index curves of risk model. (D) Uni‐Cox analyses of clinicopathologic factors and risk score with OS. (E) Multi‐Cox analyses of clinicopathologic factors and risk score with OS

Construction of nomogram

Based on risk score, age, and clinicopathological factors, a nomogram was developed for predicting the 1‐, 3‐, 5‐, and 10‐year OS incidences (Figure 4A). The calibration plots were applied to testify that whether the nomogram had an excellent concordance with the prediction (Figure 4B), which exhibited the good consistency with the actual observation.

FIGURE 4

Nomogram and calibration curves of the model. (A) Nomogram for predicting overall survival. (B) The calibration curves for 1‐, 3‐, 5‐, and 10‐year OS

The PCA and biological pathways analyses

The three‐dimensional scatter diagram of the PCA respectively showed the distribution of different patterns. The samples grouped by risk score had distinct aggregation feature (Figure 5A‐C). The results of Gene Ontology (GO) analysis include positive regulation of activation of immune response, humoral immune response, and B‐cell receptor signaling pathway (Figure 5D‐E). The results from the GSEA analysis showed different biological functions between the low‐ and high‐risk groups, such as cell cycle, DNA replication, pyrimidine metabolism, RNA degradation, spliceosome, and immune network (Figure 5F‐G). Therefore, we tried to make an immune‐related analysis based on the risk model.

FIGURE 5

The PCA and functional analyses of patients in two groups. (A‐C) The PCA 3D scatterplot of sample distribution based on necroptosis‐related lncRNAs in risk model, necroptosis‐related lncRNAs, and necroptosis‐related genes, respectively. (D‐E) The GO analysis. (F‐G) The GSEA analysis

Investigation of immune signature in risk groups

Significant differences in the immune cell infiltration were observed between the two groups (Figure 6A), and the intricate correlations existed between TIICs and 7 necroptosis‐related lncRNAs (Figure 6B). As shown in the scatter diagrams, dendritic cells activated, M0 macrophages, and M2 macrophages were positively correlated with the aforesaid risk scores, by contrast, the other TIICs had negative correlations. (Figure 6C). As for the TME score, high‐risk patients showed lower stromal scores, immune scores, and ESTIMATE score than low‐risk patients (Figure 6D).

FIGURE 6

Immune signature in two groups. (A) Expression of immune cells in the low‐ and high‐risk groups. (B) Correlations between the TIICs and 7 necroptosis‐related lncRNAs in the proposed model. (C) Correlations between risk score and immune cell types. (D) TME score in the in two groups

TMB, TIDE, and therapeutic drug sensitivity

Then, we analyzed the variations of the somatic mutations in two risk groups. The ten highest mutated genes were PIK3CA, TP53, TTN, CDH1, GATA3, MUC16, MAP3K1, MUC4, KMT2C, and PTEN. Patients in high‐risk group had markedly higher frequencies of TP53 mutation, and the opposite result was discovered with the alternation of PIK3CA and CDH1(Figure 7A–B). Compared with low‐risk group, patients in high‐risk group had higher TMB (Figure 7C). Besides, patients in the high‐risk and high‐TMB group had worst prognosis compared with the other group (Figure 7D‐E). The TIDE score was significantly lower in high‐risk group compared with low‐risk group (Figure 7F). Through drug sensitivity comparison, we found that A.443654, an AKT inhibitor, showed different IC50 between two groups, and BC patients in high‐risk group were more sensitive to this drug (Figure 7G).

FIGURE 7

TMB, TIDE, and Chemotherapeutic Sensitivity. (A‐B) The waterfall plot of somatic mutation features in two groups. (C) TMB between low‐ and high‐risk groups. (D) K–M survival curves between the high‐TMB and low‐TMB groups. (E) K‐M survival curves between the two groups. (F) TIDE score between two groups. (G) IC50 difference in A.443654

DISCUSSION

Necroptosis involved in immune responses and tumor microenvironment, and the benefits of activation of necroptosis pathways combined with immune checkpoint blockade have been demonstrated in recent study. Nowadays, it has been proved that lncRNAs engage in cancer‐related cellular pathways and have good predictive power in prognosis and diagnosis. , Emerging studies have tried to establish novel and effective lncRNAs pattern risk models in malignancy, , as well as discover molecular character and potential therapy targets of breast cancer patients. , Until now, the patterns of necroptosis‐related lncRNAs in BC and the potential capability of predicting the prognosis has not been elucidated.

In this study, we established a risk model with 7 necroptosis‐related lncRNAs, including SH3BP5‐AS1, AC012073.1, AC120114.1, LINC00377, AL133467.1, AC036108.3, and AC020663.2. Then, the patients were divided into low‐ and high‐risk groups according to the median values. The ROC and C‐index curve fatherly validated the prognostic precision of the risk score. We found that the risk score could be a yardstick for predicting prognosis. Subsequently, a nomogram was constructed for predicting prognosis of BC patients, which had an excellent concordance with the prediction.

The 3D scatterplot of the PCA showed that patients categorized by necroptosis‐related lncRNAs exhibited distinct inherent biological feature. The results of Gene Ontology (GO) analysis demonstrated that activation of immune response, humoral immune response, and B‐cell receptor signaling pathway played an important role in biological pathways. Besides, the results of GSEA showed different enrichment of genes in KEGG, including cell cycle, DNA replication, pyrimidine metabolism, RNA degradation, spliceosome, and immune network. Although breast cancer is not regarded as a highly immunogenic cancer in the past, but tumor immune microenvironment impacts on a subset of breast cancers and partial patients might be suitable to immune checkpoint blockade treatment strategies after evaluation. , Based on the above finding, we tried to make an immunity analysis by CIBERSORT and ESTIMATE method in the risk model. The risk score was positively correlated with dendritic cells activated, M0 macrophages, and M2 macrophages. ESTIMATE is a method to assess the immune cells infiltration and tumor microenvironment according to gene expression. In this study, immune scores, stromal scores, and estimate scores of high‐risk groups were significantly lower.

The somatic mutation analysis showed that the samples in high‐risk group had more frequent TP53 mutation mutated and TP53 mutations may boost immunotherapy activity in BC according to previous study. Although immunotherapy has been applied successfully in some tumor types, not all the breast cancer patients can benefit from this treatment. Therefore, it is indispensable to select appropriate biomarkers to decern the patients who are more sensitive to immunotherapy.

Hypermutated breast cancer patients may benefit from PD‐1 inhibitors, and high tumor mutation burden (TMB) is related to better therapeutic effect of immune checkpoint blockade (ICB). , In this study, the patients in the high‐risk group showed a higher TMB. Tumor immune dysfunction and exclusion (TIDE) algorithm is a method for predicting ICB response in cancer. A higher TIDE score is associated with worse ICB response and has high accuracy in predicting the survival outcome of cancer patients treated with ICB. Some recent research supported its application in predicting the therapeutic effect of ICB. , In our study, the TIDE score was significantly lower in high‐risk group. In conclusion, based on the evaluation of TMB and TIDE score, the patients in high‐risk group showed better ICB response. In addition, we used pRRophetic R package to calculate the IC50 values of chemotherapeutic drugs, and the patients with high‐risk score were more sensitive to A‐443654.

There are still limitation and weaknesses in our study. Firstly, the analysis results were not validated in vitro and in vivo, and the biological function needs to be furtherly elucidated. Secondly, in view of complexity, we did not clarify the relationship between necroptosis‐related lncRNA and tumor‐infiltrating immune cells. Thirdly, in the retrospective study, there may be some biases in the case inclusion and data processing. The collection of clinical samples and external validation will be implemented in the future.

CONCLUSION

A well‐validated risk model was constructed based on 7 necroptosis‐related lncRNAs, including SH3BP5‐AS1, AC012073.1, AC120114.1, LINC00377, AL133467.1, AC036108.3, and AC020663.2. The BC patients in the high‐risk group had worse clinical outcomes. Besides, the high‐risk patients demonstrated higher TMB and lower TIDE score, indicating the better immune checkpoint blockade response. These findings pointed novel ways of BC prognosis estimation and optimal treatment strategy.