Abstract

Background. Colorectal cancer is one of the most common gastrointestinal malignancies globally. Necroptosis has been proved to play a role in the occurrence and development of the tumor, which makes it a new target for molecular therapy. However, the role of necroptosis in colorectal cancer remains unknown yet. Our study aims to build a prognostic signature of necroptosis-related lncRNAs (nrlncRNAs) to predict the outcomes of patients with colorectal cancer and facilitate in anticancer therapy. Method. We obtained RNA-seq and clinical data of colorectal adenocarcinoma from the TCGA database and got prognosis-related nrlncRNAs by univariate regression analysis. Then, we carried out the LASSO regression and multivariate regression analysis to build the prognostic signature, whose predictive ability was tested by the Kaplan-Meier as well as ROC curves and verified by the internal cohort. Moreover, we divided the cohort into 2 groups based on median of risk scores: high- and low-risk groups. By analyzing the difference in the tumor microenvironment, microsatellite instability, and tumor mutation burden between the two groups, we explored the potential chemotherapy and immunotherapy drugs. Results. We screened out 9 nrlncRNAs and built a prognostic signature based on them. With its good prognostic ability, the risk scores can act as an independent prognostic factor for patients with colorectal cancer. The overall survival rate of patients in high-risk group was significantly higher than the low-risk one. Furthermore, risk scores can also give us hints about the tumor microenvironment and facilitate in predicting the response to the CTLA-4 blocker treatment and other chemotherapeutic agents with potential efficacy such as cisplatin and staurosporine. Conclusions. In conclusion, our prognostic signature of necroptosis-related lncRNAs can facilitate in predicting the prognosis and response to the anticancer therapy of colorectal cancer patients.

1. Introduction

Colorectal cancer is the third most common cancer in the world [1]. There are about 900,000 deaths each year, which make it the second leading cause of cancer death worldwide [2]. Due to the lack of biomarkers for early screening and prognostic prediction, many patients are not diagnosed until to an advanced stage when they cannot be surgically treated [3]. For those patients, systemic treatments are the only choice. Although the advances in diagnosis and treatment, especially the application of molecular therapeutic targeted drugs and immune checkpoint inhibitors, have relatively extended the overall survival time of patients with advanced cancer, the overall prospect for those patients is not bright. Since immunotherapy is only effective in patients with specific genotypes [4]. Therefore, it is important to identify more novel biomarkers for diagnosis and targets for anticancer therapy.

Resistance to apoptosis is one of the main hallmarks of cancer and has long been a major impediment to anticancer therapy [5, 6]. Therefore, inducing programmed death by other mechanisms is recognized as an alternative approach to overcome this obstacle [7]. Necroptosis is a caspase-independent regulatory cell death mode [8]. It could activate the mixed lineage kinase domain-like protein (MLKL) by phosphorylation signaling pathway, which is mediated by receptor-interacting protein 1/3 (RIPK1/RIPK3) [9]. Recent studies suggest that tumor cells with resistance to apoptosis may be sensitive to the necroptosis pathway [10, 11], which is expected to be a new target for anticancer therapy. As a coalescence of necrosis and apoptosis, necroptosis is regarded to play a dual role in tumors: on one hand, necroptosis can induce “necrotic” death of tumor cells, bypassing the apoptotic pathway and elicit strong adaptive immune responses via damage-associated molecular patterns (DAMPs) to stem tumor development [12]; on the other hand, necroptosis can promote tumor metastasis and progression through destruction of endothelial cells or inflammatory response [13, 14]. Additionally, necrosis-associated inflammation enhances the immunogenicity of cancer cells, which can promisingly synergize with ICBs to create new immunotherapeutic [15].

Long non-coding RNAs (lncRNAs), defined as non-protein coding RNA transcripts of over 200 nucleotides, are engaged in various cellular biological processes, including tumor progression and immune cell infiltration [16, 17]. lncRNAs play a significant regulatory part in the development of colorectal cancer. For example, some studies showed that LINC01021 [18] could affect cell cycle, proliferation, apoptosis, epithelial-mesenchymal transformation, and other processes in colorectal cancer cells, and lncRNA-p21 [19] could inhibit the invasion and metastasis of colorectal cancer cells. Moreover, lncRNAs have also been reported to be involved in tumor cell necroptosis. Tran et al. [20] found that liver cancer cells could regulate microRNA and its target genes to induce necroptosis by expressing LINC00176. However, there are few literature reports on the role of lncRNAs in the necroptosis pathway of colorectal cancer; thus, further studies are needed.

Given that, it is necessary to shed light on interactions between necroptosis-related lncRNAs and the clinicopathological characteristics, tumor microenvironment, anticancer therapy, and tumor mutation in colorectal cancer. In summary, our work can help fill the research gap on necroptosis-related lncRNAs of colorectal cancer and provide new insights into the possible pathogenesis of colorectal cancer.

2. Method and Materials

2.1. Data Collection and Preprocessing

The RNA transcriptome data (HTSEQ-FPKM format) of colorectal adenocarcinoma were obtained from TCGA-COAD and TCGA-READ projects of the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) [21], including mRNA and lncRNA expression levels of tumor samples ( =568) and normal samples ( =44), as well as corresponding clinical data (XML format), such as survival information, sex, age, TNM stage, and tumor stage. We processed the data using the limma package of R language software (version 4.1.2) and Perl language (version 5.30.3) and obtained colorectal adenocarcinoma gene expression and clinical information matrix. Simple nucleotide variation data of colorectal adenocarcinoma were also downloaded from the database, and the data type was Masked Somatic Mutation, calculated by VARSCAN software (MAF format). To control the bias, patients with overall survival (OS) of less than 30 days were excluded, and 507 colorectal adenocarcinoma patients were included. We randomly divided them into the train set and test set by R caret package, with a ratio of 1 : 1. The clinical characteristics of train set, test set, and entire set are shown as Table 1.

2.2. Selection of Necroptosis-Related Genes and lncRNAs

We obtained 67 necroptosis-related genes (NRG) from the Gene Set Enrichment Analysis (GSEA) database (http://www.gsea-msigdb.orggseaindex.jsp) [22, 23] and previous literature [24]. Pearson’s correlation analysis was performed (|Pearson R|>0.5, ) to get necroptosis-related lncRNAsm (nrlncRNAs), using R limma package. We furtherly conducted the differential analysis (Log2fold change (FC)>1, false discovery rate (FDR)<0.05, and ) to identify the nrIncRNAs that were significantly differentially expressed between tumor and normal samples.

2.3. Establishment and Validation of the Prognostic Signature

First, the necroptosis-related lncRNAs that are closely related to prognosis were screened out by univariate Cox regression analysis (). Then, to avoid overfitting and ensure the minimum amount of lncRNAs as well as complete information, the least absolute shrinkage and selection operator (LASSO) regression was performed with 10-fold cross-validation and 1000 random circles. Finally, we conducted multivariate Cox regression to identify the nrlncRNAs that could serve as independent prognostic factors for modeling. The following formula was devised based on the risk model to calculate the risk score of each patient: where the n was the number of lncRNAs, was the expression level of lncRNAs, and was the correlation coefficient between lncRNAs and survival data. Based on the median risk score of the train set, subgroups were established including low-risk and high-risk groups. The network between nrlncRNAs and NRGs was drawn using Cytoscape software (V3.8.0) and the survival curve and receiver operator characteristic (ROC) curve were established through R survival, survminer, and timeROC package to evaluate the predictive performance of the signature.

2.4. Clinical Applicability of the Signature and Construction of Nomogram

We drew the survival curves to compare the differences in survival between different clinical subgroups based on previously obtained clinical information. The risk scores were verified by univariate and multivariate COX regression to decide whether it is an independent prognostic factor. ROC curves were utilized to compare the predictive efficacy of risk scores and different clinical variables. Furthermore, the nomogram was made, using R rms package to predict the 1-, 3-, and 5-year survival rate, with age, gender, TNM stage, and risk score as variables. Meanwhile, the Hosmer-Lemeshow test was carried out to test the consistency of the predicted results with the actual situation, and decision curve analysis (DCA) was performed to compare the applicability of different clinical variables and nomograms.

2.5. Gene Set Enrichment Analyses

We conducted the gene enrichment analysis to uncover the biological processes that might differ significantly between the high- and low-risk groups. The Gene Ontology (GO) Gene Set (c5.go.v7.5.1.symbols.gmt) was downloaded from the Gene Set Enrichment Analysis website (http://www.gsea-msigdb.org/), including cellular component (CC), molecular function (MF), biological process (BP) gene set data, and Kyoto Encyclopedia of Genes and Genomes (KEGG) (KEGG.v7.4. symbols. GMT) gene set data. We analyzed the enrichment of two risk groups in different pathways and functions using GSEA software (version 4.1.0), whose differences were statistically significant when and FDR <0.25.

2.6. Tumor Microenvironment and Immune Checkpoints

In order to know more about the tumor microenvironment (TME), we calculated the stromal score, immune score, and ESTIMATE score of the samples to deduce the tumor purity via the R estimate package. We then downloaded the immune cell infiltration files from the Timer 2.0 [25] database (http://timer.cistrome.org/). The results of the immunization assessment calculated by Timer, CIBERSORT, XCELL, QUANTISEQ, MCPcounter, and EPIC software platforms were included, and a correlation analysis was conducted between the results and risk score (). Finally, we evaluated the distribution of immune cells and immune function in two groups and compared the expression differences of immune checkpoints between the high- and low-risk groups by performing single sample Gene Set Enrichment Analysis (ssGSEA) and -test. The above results were presented by box plot, bubble plot, and relationship diagram using R ggplot2, ggpubr, and packages.

2.7. Immunotherapy and Chemotherapeutics Sensitivity Analysis

To evaluate the performance of the prognostic signature in guiding immunotherapy, we downloaded the immunophenoscore (IPS) data of colorectal adenocarcinoma from the Cancer Immune Altas (TCIA, https://tcia.at/home) [26] database to compare the responses of patients to PD-1 and CTLA-4 blockers treatment in high- and low-risk groups. Moreover, we compared the half-maximal inhibitory concentration (IC50) of different chemotherapy drugs to identify the promising therapeutic substances, of which data were obtained from Genomics of Drug Sensitivity in Cancer database (GDSC) (https://www.cancerrxgene.org/) and analyzed by oncoPredict package [27].

2.8. Microsatellite Instability and Tumor Mutation Burden

Similarly, we obtained the data related to microsatellite instability (MSI) in patients with colorectal adenocarcinoma from the TCIA database and used R ggplot package to draw the percentage chart of microsatellite stability (MSS), MSI-L, and MSI-H in high- and low-risk groups. A waterfall map was drawn through R maftools package to display the mutation frequency and mutation type of tumor mutation genes in the two risk groups based on the previously downloaded data. Ultimately, we compared the tumor mutation burden (TMB) and its effect on the prognosis by difference analysis and survival curve. TMB is determined by the ratio of the frequency of somatic mutations to the length of exon effective regions [28], and the somatic mutation data were derived from TCGA-COREAD mutation data.

3. Results

3.1. Identification of Necroptosis-Related lncRNAs

The entire analysis process is shown in Figure 1. We obtained 568 colorectal adenocarcinoma samples and 44 normal samples from the TCGA database and sorted out the expression matrix of 56461 lncRNAs and 67 necroptotic-related genes (NRG) corresponding to each sample. 1425 necroptosis-related IncRNAs (nrlncRNAs) were identified by Pearson’s correlation analysis (|Pearson R|>0.5, ). The Sankey diagram (Figure 2(a)) showed the connections between 67 NRGs and nrlncRNA gene sets. Through differential analysis (LogFC>1, DFS<0.05), we found 747 significantly differentially expressed nrlncRNAs between tumor and normal samples. The volcano plot (Figure 2(b)) showed that 679 nrlncRNAs were significantly up-regulated and 68 nrlncRNAs significantly down-regulated in tumor samples.

3.2. Establishment and Verification of Prognostic Signature

In order to screen out nrlncRNAs that can predict the prognosis of patients, we used univariate COX regression analysis to obtain 19 nrlncRNAs that significantly affected the overall survival rate (OS) of patients (all ). The forest map (Figure 3(a)) showed the hazard ratio (95% confidence interval) of these 19 nrlncRNAs, and the heat map (Figure 3(b)) reflected their expression levels in tumor or normal samples. To reduce the overfitting combination of linear regression and improve the model’s accuracy, we adopted LASSO regression to reduce dimensionality and got the variables corresponding to the minimum partial likelihood deviance (Figure 3(c)), namely, the identified 13 nrlncRNAs. Through multivariate COX regression, 9 nrlncRNAs (Figure 3(g)) that could be used as independent prognostic factors were finally obtained for modeling. The Sankey diagram (Figure 3(e)) and network diagram (Figure 3(f)) showed the corresponding connections between these nrlncRNAs and necroptosis mRNAs. The risk scores based on the expression level of nrlncRNAs was as follows: Risk Score = MYOSLID expi× (3.2358) + AC006111.2 expi× (-1.9913) + AC245100.5 expi × (0.5341) + AL161729.4 expi × (1.1819) + AL355312.2 expi × (-3.0215) + AL137782.1 expi × (-0.8773) + NSMCE1-DT expi × (2.9534) + LINC02257 expi × (0.8564) + LINC00513 expi × (-0.3655).

The median risk score of patients in the train set acted as the node to divide the cohort into high- and low-risk groups, while the test set was used for in-cohort verification. Figure 4 shows the nrlncRNAs expression level, risk score, and survival time distribution of high- and low-risk groups in the entire set, train set, and test set. Survival curve indicated that overall survival rate in the high-risk group was significantly lower than low-risk group (entire set: , train set: , test set: ). The ROC curve showed that the one-year prognosis area under the curve (AUC) was 0.751, 0.860, and 0.657, three years 0.767, 0.893, and 0.646, and five years 0.751, 0.860, and 0.657 in the three cohort sets, respectively. In addition, the cohort was divided into 9 pairs of high- and low-expression groups based on each of the selected nrlncRNAs, respectively, and 7 nrlncRNAs showed significant differences in overall survival rate (Supplementary 1).

3.3. Clinical Applicability of Prognostic Signature

It was found that risk scores could serve as an independent risk factor for prognosis among clinical factors including age, gender, tumor stage, and TNM stage (Figures 5(a) and 5(b)), by univariate COX regression analysis and multivariate COX regression analysis, with respective hazard ratio 1.055 (95% CI: 1.039-1.072) and 1.051 (95% CI: 1.034-1.069). ROC curve (Figure 5(c)) showed that the area under the signature curve was 0.751, greater than the 0.741 of the clinical staging, which was widely used as a prognostic predictor in clinical practice. According to the heat map (Figure 5(d)), the expression levels of nrlncRNAs for the signature varied greatly in different tumor stages or TNM stages. Furthermore, we compared the Kaplan-Meier curves of the high- and low-risk groups under different clinical conditions (Figure 5(e)) and found that the overall survival rate of the high-risk group was far worse than that of the low-risk group, except in the case of early T stage (T1-T2) and distant metastasis (M1), possibly due to the insufficient number of patients. All this indicated that the prognostic signature applied to patients with different clinical conditions.

3.4. Construction and Verification of Nomogram

Using variables such as age, sex, clinical stage, TNM stage, and risk score, we constructed a nomogram that predicted 1-, 3-, and 5-year survival. In Figure 6(a), a 67-year-old female CRC patient with stage T3N0M0 II was grouped into the high-risk group with a total score of 197. The predicted 1-year survival rate of this patient was 88.7%, 3-year survival rate 71.2%, and 5-year survival rate 49.8%. Principal component analysis (PCA) diagram (Figure 6(b)) showed that it is easy to distinguish the patients between high- and low-risk groups. We also used 1-, 3-, and 5-year calibration curves (Figure 6(c)) to demonstrate that the nomograms were broadly consistent with the observed OS. Decision curve analysis (Figure 6(d)) showed that the nomogram was more effective than other clinical factors in predicting patient outcomes.

3.4.1. Gene Set Enrichment Analysis

Compared with the low-risk group, genes in the high-risk group were significantly enriched in multiple biological functions and pathways, including extracellular matrix structure and some inflammatory pathways (Figure 6(e)), which justified the necessity of further research in the tumor microenvironment and associated immunity. However, no significant gene enrichment was found in the low-risk group. Additionally, we found that some high-risk genes were significantly enriched in tumor pathways and immune pathways, such as MAPK, Wnt, TGF-β, Toll-like receptor signaling pathway, and the like.

3.5. Tumor Microenvironment and Immune Infiltration

Using the ESTIMATE algorithm, we found that the stromal score, immune score, and estimate score in the high-risk group were significantly higher than those in the low-risk group (Figure 7(a)). In terms of immune cell infiltration (Figure 7(b)), the expression data of immune cells on different software platforms showed that the number of most immune cells infiltrated was positively correlated with risk scores, including cancer-associated fibroblast (EPIC), Macrophage M2 (CIBERSORT_ABS), T Cell Regulatory (Tregs) (QUANTISEQ), Myeloid Dendritic Cell (TIMER), and the like, while the number of some immune cells infiltrated was negatively correlated with the risk scores, such as NK cell (QUANTISEQ) and T cell CD4+ Th1 (XCELL). Similarly, ssGSEA scores (Figure 7(c)) of multiple immune functions were significantly higher in the high-risk group than in the low-risk group.

3.6. Immune Checkpoints and Anticancer Therapy

The expression levels of 24 immune checkpoints were significantly higher in the high-risk group than the low-risk group (Figure 7(d)), including CD274 (PD-1) and other novel immune checkpoints such as NRP1, BTLA, TIGIT, and VTCN1. It indicated that the high-risk group might have greater immunotherapy promises. As shown in Figure 8(b), the half-maximal inhibitory concentration (IC50) of therapeutic substances in the high-risk group, including cisplatin, staurosporine, and cyclophosphamide, was notably lower than that of the low-risk group, which showed that patients in the high-risk group were more sensitive to these compounds. However, there were no significant differences in drug sensitivity to 5-fluorouracil, oxaliplatin, and irinotecan between the two groups. Furthermore, by comparing the immunophenoscore (IPS) of two groups when treated with PD-1 blocker and CTLA-4, we found that the low-risk group had a better response to the latter, while no significant difference was observed in response to PD-1 blocker or the combination of the two between the two groups.

3.7. Microsatellite Instability and Tumor Mutation

Compared with the low-risk group, the high-risk group had larger proportions of microsatellite instability-high (MSI-H), microsatellite instability-low (MSI-L) (19% vs.12%, 18% vs.15%), and lower microsatellite stability (MSS) (63% vs.73%) (Figure 7(c)), while there was no significant difference in risk scores between MSS/MSI-L and MSI-H patients. In the two risk groups, APC, TP53, TTN, KRAS, and SYNE1 were the top five mutation genes, of which the proportion of APC was higher in the low-risk group, and the proportion of TP53, TTN, KRAS, and SYNE1 in the high-risk group was higher (Figure 7(d)). The high-risk group’s tumor mutation burden (TMB) was higher than that of the low-risk group, although there was no significant statistical difference. The TMB remarkably impacted the overall survival rate. Survival was worse in the TMB-H group than TMB-L group at 1 to 5 years, but at 5 to 7 years, it was the other way around. However, the effect appeared weaker than risk scores, for patients in the high-risk group had much worse outcomes than those in the low-risk group, even at a low tumor low mutation burden.

4. Discussion

Colorectal cancer is a highly malignant digestive system tumor with high morbidity and mortality. Although a systemic treatment protocol has been formed on the basis of surgery, chemotherapy, and radiotherapy, the prognosis is still not optimistic. One of the reasons is chemotherapy resistance, directly or indirectly [29]. Necroptosis is considered to be an alternative for apoptosis, which plays a role in inducing necrotic cell death when tumor cells develop drug resistance due to dysregulation of the apoptosis mechanism [30]. Moreover, more and more studies have shown that necroptosis also got involved in tumorigenesis.

Given that, Huang et al. [31] used necroptosis-related miRNAs to predict the prognosis of colon cancer patients, which have shown fair performance. However, we only have a limited understanding of the interactions between lncRNAs and necroptosis in colorectal cancer, which requires further study. This study identified 9 necroptosis-related lncRNAs (nrlncRNAs) from the public database, established a prognostic signature, and grouped them by risk scores to explore the signature’s performance in predicting immune infiltration and guiding immunotherapy. Our study suggested that prognostic signature based on necroptosis-related lncRNAs can be used for prognostic stratification of colorectal cancer patients and help to decipher the molecular mechanism of colorectal cancer.

Previous researches have proved that necroptosis-related lncRNAs are closely related to gastric cancer patient’s prognosis and can be used to distinguish between cold and hot tumors and guide immunotherapy [24]. This has demonstrated nrlncRNAs’ research potential in gastrointestinal tumors. Through statistical analysis and screening, we managed to obtain 9 nrlncRNAs that can be utilized to predict prognosis, including MYOSLID, AC006111.2, AC245100.5, AL161729.4, AL355312.2, AL137782.1, NSMCE1-DT, LINC02257, and LINC00513, of which the former five lncRNAs were risk factors and the latter four were protective factors. MYOSLID, with the highest hazard ratio, has the most significant correlation with survival rate (coefficient =3.236, hazard ratio =16.54 (2.92-93.64), ), which has been reported to regulate the biological processes of various tumor cells. Han et al. [32] demonstrated the critical role of the MYOSLID-miR-29c-MCL-1 axis in the tumorigenesis of gastric cancers. Xiong et al. [33] found that MYOSLI can promote the invasion and metastasis of squamous cell carcinomas of the head and neck by regulating epithelial-mesenchymal transformation. Moreover, LINC02257 [34] has also been used for predicting the outcomes of colorectal cancers, whose expression level is associated with the immune state. LINC00513 [35] is considered to play a part in systemic lupus erythematosus and acts as a positive regulator of interferon signaling pathways. Other nrlncRNAs are not reported in relevant literature, whose biological effects need to be explored. There are close interactions between these nrlncRNAs and necroptotic genes, indicating that they share some common mechanisms in necroptotic pathways. Although the expressions of these lncRNAs are all significantly positively correlated with necroptotic genes, they have different effects on prognosis, which is consistent with the contradiction of necroptosis in tumors.

The prognostic signature based on these nrlncRNAs has good predictive performance. Clinically, the most widely used prognostic indicators of colorectal cancer are TNM stage and clinical stage, while our signature is superior to these indicators and applicable to different ages, genders, and stages, whose effectiveness is furtherly verified by the internal cohort. Besides, we also drew a nomogram that included risk scores and clinical indicators to facilitate the prediction of the 1-, 3-, and 5-year survival rate of a single patient, whose prediction results had good consistency with the actual situation.

Gene enrichment analysis revealed that genes in the high-risk group were significantly enriched in common signaling pathways of tumors and immune signaling pathways, which may be the underlying molecular mechanisms behind the worse prognosis of patients in the high-risk group. In addition, genes in the high-risk group are also involved in extracellular matrix related biological processes. Extracellular matrix is a protein network surrounding normal cells and cancer cells, acting as an important component of tumor microenvironment, which greatly influences tumor proliferation, angiogenesis, adhesion, movement, invasion, and metastasis [36, 37]. Thus, it is necessary to probe into the tumor microenvironment.

Immune cell infiltrated in tumor tissue could influence the onset and development of tumor and its response to immunotherapy [3841]. Overall, the high-risk group is burdened with higher tumor cell purity and more immune cell infiltration than the low-risk group. Specifically, the infiltration of M2-type macrophages, Tregs, and dendritic cells is significantly positively correlated with the risk scores. Studies have shown that massive M2-type macrophages in tumor tissues signified poor prognosis [42]. On one hand, Tregs cells exert greater immunosuppressive function by upregulating the expression of cytotoxic T lymphocyte antigen-4 (CTLA-4) and programmed death receptor-1 (PD-1) in colorectal cancer [43]. On the other hand, the migration of T cells to the tumor is inhibited by Treg cells by regulating chemokines [44] and secreting cytokines that block antitumor immunity [45]. The degree of dendritic cell infiltration is positively correlated with the metastasis and stage of colorectal cancer [46, 47]. In part, NK cells and Th1 cells are negatively correlated with the risk scores, which perform an antitumor function [48, 49].

Immune checkpoint inhibitors (ICIs) enhance the cytotoxicity effect of T cells on tumor cells by acting on co-inhibitory receptors such as CTLA-4 and PD-1 or their ligands such as programmed death ligand-1 (PD-L1) [50]. The application of ICIs in colorectal cancer is beneficial in a minority of patients with high immunogenicity. In our study, the low-risk group had better response to the CTLA-4 blocker treatment, even though there is a larger proportion of MSI-H patients. In the high-risk group, more immune checkpoint genes were expressed. The high exposure of these immune checkpoints and infiltration of M2 macrophages, Treg cells, and MODS cells constitute the immunosuppressive environment in the high-risk group. Research has shown that necroptosis-induced CXCL1 and Mincle signaling promote macrophage-induced adaptive immune suppression enabling pancreas oncogenesis [51].

Therefore, the five lncRNAs highly expressed in the high-risk group may be involved in the regulation of necroptosis-induced immunosuppression in colorectal cancer, which may account for the patients’ relative insensitivity to ICIs in the high-risk group.

Recent studies have found that some compounds can play a role in anticancer therapies by inducing necroptosis, which may provide us a new promising therapeutic approach for bypassing the acquired or intrinsic apoptosis resistance. For example, cisplatin can induce necroptosis in esophageal and lung cancer cells [52, 53] without being affected by apoptosis resistance. Staurosporine has been reported to induce RIPK1 and MLKL-dependent necroptotic cell death in leukemia cells when caspase activation is compromised [54]. In this study, the patients in the high-risk group were more sensitive to cisplatin and staurosporine, which implies the value of both compounds in the treatment of colorectal cancer and the role of lncRNAs in the regulation of necroptosis induction. However, the guidance for the use of chemotherapy drugs including 5-fluorouracil, oxaliplatin, irinotecan, and raltitrexed is limited, which are commonly used in colorectal cancer.

Tumor mutations are not rare in colorectal cancers, and the tumor cells with a high tumor mutation burden (TMB) express more tumor antigens, which intensifies the immune-killing effect [55]. However, TMB did not differ significantly between the high- and low-risk group, suggesting similar immunogenicity due to mutation exposure between the two groups. Moreover, the tumor mutation burden also has an impact on prognosis. Studies [56] have shown that patients with high TMB have a longer overall survival time compared with colorectal cancer patients with low TMB, which is consistent with the long-term survival prospects in our cohort.

Our study had some limitations: Firstly, the prognostic signature and the nomogram lacked external queue verification [57]. We tried to download gene chips from the Gene Expression Omnibus (GEO) database for external verification, but the expression level of key lncRNAs could not be extracted due to the incomplete sequencing genes. Secondly, only limited data analysis cannot fully elucidate the specific regulatory role of 9 nrlncRNAs in necroptosis of colorectal cancer and further exploration in vivo and in vitro experiments is needed. Finally, our data mainly came from the TCGA database, and the missing data is not random, which may cause some bias.

5. Conclusion

In conclusion, our study established a relatively reliable prognostic signature on the basis of necroptosis-related lncRNAs, which facilitates in predicting the prognosis as well as tumor immune microenvironment, and guiding the antitumor treatment of colorectal cancer patients according to the risk stratification. Our study also showed that lncRNAs may get involved in the double nature of necroptosis in colorectal cancer, and its molecular mechanisms require further experiment verification.

Abbreviations

AUC:Area under the curve
BP:Biological process
CC:Cellular component
CRC:Colorectal cancer
CTLA-4:Cytotoxic T lymphocyte antigen-4
DCA:Decision curve analysis
GEO:Gene Expression Omnibus
GO:Gene Ontology
GSEA:Gene Set Enrichment Analysis
GDSC:Genomics of Drug Sensitivity in Cancer database
IPS:Immunophenoscore
KEGG:Kyoto Encyclopedia of Genes and Genomes
LASSO:Least absolute shrinkage and selection operator
LncRNA:Long non-coding RNA
MSI:Microsatellite instability
MSS:Microsatellite stability
MLKL:Mixed lineage kinase domain-like protein
MF:Molecular function
NRG:Necroptosis-related gene
nrlncRNA:Necroptosis-related lncRNA
OS:Overall survival
PCA:Principal component analysis
PD-L1:Programmed death ligand -1
PD-1:Programmed death receptor-1
ROC:Receiver operator characteristic
RIPK1:Receptor-interacting protein 1
RIPK3:Receptor-interacting protein 3
ssGSEA:Single sample Gene Set Enrichment Analysis
TCGA:The Cancer Genome Atlas
TCIA:The Cancer Immune Altas
IC50:The half-maximal inhibitory concentration
TME:Tumor microenvironment
TMB:Tumor mutation burden.

Data Availability

Data is available at TCGA database (https://portal.gdc.cancer.gov/).

Disclosure

The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

The author contributions were presented as follows: (I) Hanyu Xiao contributed to the conception and design; (II) Hong Chen contributed to the administrative support; (III) Yong Wang contributed to the provision of study materials or patients; (IV) Suhe Lai contributed to the collection and assembly of data; (V) Qidan Pang contributed to the data analysis and interpretation; (VI) all authors contributed to the manuscript writing; (VII) all authors contributed to the final approval of the manuscript. Hanyu Xiao and Qidan Pang contributed equally to this work and shared the first authorship.

Supplementary Materials

Supplementary 1. The survival curve of patients with high-expression and low-expression of 9 identified necroptosis-related lncRNAs.

Supplementary 2. Potential drugs prediction for high-risk group versus low-risk group.