Abstract Necroptosis is intimately associated with the initiation and progression of colon adenocarcinoma (COAD). However, studies on necroptosis-related genes (NRGs) and the regulating long non-coding RNAs (NRGlncRNAs) in the context of COAD are limited. We retrieved the cancer genome atlas (TCGA) to collect datasets of NRGs and NRGlncRNAs on COAD patients. The risk model constructed using Cox and least absolute shrinkage and selection operator (LASSO) regression was then employed to identify NRGs and NRGlncRNAs with prognostic significance. Subsequently, we validated the results using gene expression omnibus (GEO) datasets from different populations, conducted Mendelian randomization (MR) analysis to explore the potential causal relationships between prognostic NRGs and COAD, and conducted cell experiments to verify the expression of prognostic NRGlncRNAs in COAD. Furthermore, we explored potential pathways and regulatory mechanisms of these prognostic NRGlncRNAs and NRGs in COAD through enrichment analysis, immune cell correlation analysis, tumor microenvironment analysis, immune checkpoint analysis, tumor sample clustering, and so on. We identified eight NRGlncRNAs (AC245100.5, [34]AP001619.1, LINC01614, [35]AC010463.3, [36]AL162595.1, ITGB1-DT, LINC01857, and LINC00513) used for constructing the prognostic model and nine prognostic NRGs (AXL, BACH2, CFLAR, CYLD, IPMK, MAP3K7, ATRX, BRAF, and OTULIN) with regulatory relationships with them, and their validation was performed using GEO and GWAS datasets, as well as cell experiments, which showed largely consistent results. These prognostic NRGlncRNAs and NRGs modulate various biological functions, including immune inflammatory response, oxidative stress, immune escape, telomere regulation, and cytokine response, influencing the development of COAD. Additionally, stratified analysis of the high-risk and low-risk groups based on the prognostic model revealed elevated expression of immune cells, increased expression of tumor microenvironment cells, and upregulation of immune checkpoint gene expression in the high-risk group. Finally, through cluster analysis, we identified tumor subtypes, and the results of cluster analysis were essentially consistent with the analysis between risk groups. The prognostic NGRlncRNAs and NRGs identified in our study serve as prognostic indicators and potential therapeutic targets for COAD, providing a theoretical basis for the clinical diagnosis and treatment of COAD and offering guidance for further research. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-024-73168-3. Keywords: Necroptosis, lncRNA, Colon adenocarcinoma, Prognosis, Mendelian randomization, Risk model Subject terms: Cancer, Oncology, Risk factors Introduction Colorectal cancer is among the most prevalent malignant tumors of the digestive tract globally^[37]1. In 2020, nearly 1.93 million new cases of colorectal cancer were diagnosed worldwide, with approximately 930,000 deaths. The incidence and mortality rate of colorectal cancer ranks third and second, respectively, among all types of cancer worldwide^[38]2. Colorectal cancer consists of two major categories: colon cancer and rectal cancer. The primary histological types of colon cancer are adenocarcinoma, mucinous adenocarcinoma, and undifferentiated carcinoma, with adenocarcinoma being the most common type. Colon adenocarcinoma (COAD) is a common malignancy of the colonic glandular epithelium, with its incidence rising among individuals under 50^[39]3–[40]5. while early-stage COAD can often be effectively treated with surgical intervention, improving the 5-year survival rate, most cases are diagnosed at an advanced stage. At this stage, treatment outcomes are suboptimal, with a high incidence of malignant recurrence and distant metastasis, resulting in a low 5-year survival rate^[41]6,[42]7. Therefore, it is essential to explore the pathogenesis of COAD, investigate its gene regulation patterns, identify novel prognostic markers and potential treatments, and ultimately improve patient survival. Necroptosis, a form of programmed cell death, helps the host defend against certain pathogenic invasions^[43]8,[44]9. In addition, necroptosis plays a crucial role in the initiation and progression of cancer, influencing tumorigenesis, tumor immunity, and metastasis^[45]10,[46]11. However, its impact on tumors is twofold. On the one hand, necroptosis may promote tumor metastasis and progression^[47]12–[48]14. Conversely, when apoptotic pathways are compromised, necroptosis may inhibit tumor progression^[49]15–[50]17. Studies suggest that necroptotic cells regulate tumor regrowth through the receptor-interacting serine/threonine kinase 1 (RIP1) / receptor-interacting serine/threonine kinase 3 (RIP3) / mixed lineage kinase domain-like pseudokinase (MLKL) / c-Jun N-terminal kinases (JNK) / C-X-C motif chemokine ligand 8 (IL-8) signaling pathway. MLKL/JNK/IL-8 as drug candidates may help prevent tumor regrowth, enhancing the effectiveness of radiotherapy in colorectal cancer treatment^[51]18. An in-depth understanding of the regulatory mechanisms between necroptosis and COAD has significant implications for developing new anti-cancer therapies. Long non-coding ribonucleic acids (lncRNAs) are RNA transcripts longer than 200 nucleotides that do not have protein-coding functions^[52]19. LncRNAs are involved in various disease processes, including regulating tumor development, metastasis, and invasion^[53]20–[54]22. Current research indicates that lncRNAs play a crucial role in the occurrence and progression of COAD. For instance, the lncRNA Pvt1 oncogene (PVT1) promotes COAD progression through the PVT1/miR-30d-5p/RUX2 competitive endogenous RNA (ceRNA) network; lncRNA LINC00662 promotes the progression and metastasis of COAD by competitively binding with miR-340-5p to regulate the co-expression of claudin 8 (CLDN8) / interleukin 22 (IL22) and the activation of the ERK signaling pathway; lncRNA H19 can serve as a biomarker aiding in the detection of COAD^[55]23–[56]25. These findings reveal the molecular mechanisms linking lncRNAs to COAD and highlight new biomarkers and therapeutic targets for the disease. However, current research on necroptosis-related genes (NRGs) and their regulatory long non-coding RNAs (NRGlncRNAs) in COAD remains insufficient and lacks highly credible validation. Therefore, this study aims to identify NRGs and NRGlncRNAs by analyzing COAD datasets from the cancer genome atlas (TCGA) database and constructing a COAD prognostic model of NRGlncRNAs to screen for prognostic NRGs and NRGlncRNAs in COAD. Subsequently, the GEO datasets and cell experiments will be used to validate the screening results, and potential causal relationships between prognostic NRGs and COAD will be analyzed using Mendelian randomization (MR) analysis. Finally, the roles of prognostic NRGlncRNAs in COAD will be further analyzed through various methodologies. This study aims to identify prognostic NRGlncRNAs and NRGs as potential indicators and therapeutic targets for COAD, providing valuable insights for clinical diagnosis and treatment. The sections of the manuscript are summarized as follows: section “[57]Materials and methods” introduces the data sources and methods of this study, focusing on constructing the prognostic risk model using Cox and least absolute shrinkage and selection operator (LASSO) regression analyses. Section “[58]Results” describes the results obtained in this study, focusing on prognostic NRGlncRNAs and NRGs. Section “[59]Discussion” provides an analysis and discussion of the findings. Section “[60]Conclusion” summarizes the study’s significance. The research design is shown in Fig. [61]1. Fig. 1. [62]Fig. 1 [63]Open in a new tab Flow chart of our study. Materials and methods Collection and processing of datasets We downloaded gene expression, transcriptome, and clinical data for COAD from TCGA ([64]https://portal.gdc.cancer.gov/) (Supplementary 1—Clinical information of TCGA clinical data) and GEO ( [65]https://www.ncbi.nlm.nih.gov/geo/) databases. Specifically, the transcriptome data from TCGA were decompressed to form a matrix of 39 normal tissues and 398 tumor tissues, and the Ensembl IDs were converted to gene symbols (The genome annotation file was obtained from [66]https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.39/). Perl code was used to extract clinical data, including age, gender, survival time, survival status, grade, stage, and TNM staging information. Next, based on the reference document “GRCh38.p13” ([67]https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39/), we distinguished between mRNA and lncRNA. Cases with survival times less than 30 days were removed to prevent potential biases. We obtained two datasets from the GEO database. The gene expression data of [68]GSE44076 are derived from tumor samples of 98 COAD individuals, paired adjacent normal mucosa, and 50 healthy colon mucosa samples. The lncRNA expression data of [69]GSE70880 are derived from tumor and adjacent normal tissue samples of 20 COAD patients. The lncRNA data were re-annotated using gene symbols. We collected 67 NRGs from various reviews and reports (Supplementary 2—Necroptosis Related Genes). Summary data for prognostic NRGs’ genome-wide association studies (GWAS) were collected from the integrated epidemiology unit (IEU) database ([70]https://gwas.mrcieu.ac.uk/) as exposure factors. The outcome data of COAD was collected from the FinnGen ([71]https://www.finngen.fi/en/access_results), which includes 218,792 participants of European ancestry (1,396 cases and 217,396 controls). These data were used to analyze the potential causal relationship between prognostic NRGs and COAD through MR. The characteristics of the datasets used in this study are detailed in Supplementary 3. Analysis of co-expression between NRGs and lncRNAs, and differential expression analysis of NRG-lncRNAs The expression levels of NRGs in each sample were extracted using the “limma” package in R. Correlation tests were then performed on the expressions of NRGs and lncRNAs in tumor samples^[72]26. Regulatory relationships were identified based on a correlation coefficient of ≥ 0.4 or ≤ -0.4, and a P value of ≤ 0.001. A co-expression network was constructed to visualize these results. The differential expression of NRGlncRNAs between tumor and normal samples was analyzed using thresholds of |logFC| ≥1 and FDR ≤ 0.05. The results were visualized using heatmaps and volcano plots generated with R packages. Construction of the NRG-lncRNA prognostic model, prognosis analysis, and visualization The datasets containing univariate significant NRGlncRNAs expression data and clinical survival data were randomly and evenly split into a train group and a test group (1:1)^[73]27. A Lasso regression model was then constructed and analyzed for NRGlncRNAs expression, survival time, and survival status. Cross-validation was conducted, and the NRGlncRNAs with the smallest errors were selected as the significant NRGlncRNAs for the Lasso regression. A Cox regression model was constructed by combining the significant NRGlncRNAs and their expression levels. The model calculation formula was obtained using the train group, where the NRGlncRNA expression levels were multiplied by their coefficient and summed. The risk scores were calculated using this formula, dividing the Train group into low-risk and high-risk categories, with the median value serving as the threshold. The median was chosen as the cut-off point because its statistical properties are less affected by outliers than the mean. This method ensures a balanced division of the cohort, which aids in assessing the prognostic utility of the risk scores. Any given patient’s tumor sample can be directly compared with the median risk score once the expression of these lncRNAs has been determined and processed through the risk scoring formula. This helps ascertain whether the patient falls into the low-risk or high-risk category, facilitating personalized prognostic assessment and potential treatment customization. The same process is applied to the test group, where risk scores are calculated based on the formula and samples are divided into high-risk and low-risk groups according to the median value from the train group^[74]28. The subsequent analysis yielded prognostic NRGlncRNAs, and differential analysis was conducted. We also performed GO (Gene ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses on these prognostic NRGlncRNAs. We further analyzed the regulatory relationships between NRGs and prognostic NRGlncRNAs, identifying prognostic NRGs with regulatory relationships to prognostic NRGlncRNAs. Patients in the train and test groups were categorized into high and low-risk groups based on their risk scores. Survival curves, risk curves, and receiver operating characteristic (ROC) curves were generated for independent prognostic analysis. Patients were also stratified into different groups based on age, gender, and stage, and the model was validated in these clinical subgroups. Finally, a nomogram was generated incorporating patients’ clinical features, and the prognosis was predicted based on the prognostic scores of different features. Validation of prognostic NRGs and NRG-lncRNAs using GEO datasets We used R packages such as “limma” and “ggpubr“([75]https://github.com/kassambara/ggpubr) to perform differential expression analysis of prognostic NRGs and NRGlncRNAs between normal and COAD samples from the GEO databases. This analysis yielded significantly differentially expressed prognostic NRGs and NRGlncRNAs. These differentially expressed genes were visualized using box plots and heatmaps, providing multi-dataset validation of these prognostic NRGs and NRGlncRNAs. Causal relationships between prognostic NRGs and COAD using MR analysis To investigate the causal relationships between prognostic NRGs and COADs, we conducted two-sample MR analyses using the “TwoSampleMR” package^[76]29. This package supports two-sample MR analysis, utilizing single-nucleotide polymorphism (SNP) data for exposures and outcomes from different GWAS, and outputting estimates of the causal effect of the exposure on the outcome, along with their statistical significance, to evaluate the causal relationship^[77]30. To maximize the number of SNPs for sufficient prognostic NRGs, we set the following parameters: (1) The P-value for the genetic instrument in association with the exposure should be less than 5 × 10^-5; (2) When conducting linkage disequilibrium (LD) clumping, we set the r² threshold to 0.4; (3) When conducting clumping analysis, we set the window size to 10,000 kilobase pairs;4) F-statistic for SNPs related to the exposure factor should be > 10; 5) We apply Bonferroni correction to adjust the significance level threshold^[78]31–[79]33. Gene set enrichment analysis (GSEA), immune cell correlation analysis, single-sample gene set enrichment analysis (ssGSEA), tumor microenvironment (TME) analysis, checkpoint analysis, and drug sensitivity screening We imported expression data and clinical data into the GSEA V.4.2.1 software (Which uses the GSEA computational method to determine whether a predefined set of genes shows statistically significant, concordant differences between two phenotypes. [80]https://www.gsea-msigdb.org/gsea/index.jsp) to perform pathway enrichment analysis, and select the five pathways with the smallest P-values in the high-risk and low-risk groups for presentation^[81]34. Next, we conducted an immune cell correlation analysis, calculating correlation coefficients between immune cells and risk scores, which were visualized as bubble plots using R. Following this, we performed ssGSEA analysis using the “GSEABase” R package (Which provides classes and methods to support GSEA, including gene set representation, gene-phenotype association, gene set management, identifier type conversion, and gene set creation. 10.18129/B9.bioc.GSEABase) (The LM22 immune gene set can be obtained from [82]https://www.nature.com/articles/nmeth.3337#MOESM207), comparing differences in immune cells and immune functions between high-risk and low-risk groups. Subsequently, the “estimate” package was used to score the tumor microenvironment of the samples, including stromal cell scores, immune cell scores, and ESTIMATE scores (the sum of the first two), and differences in scores of samples in high-risk and low-risk groups were compared^[83]35. Then, a differential analysis of immune checkpoints in high-risk and low-risk groups was conducted, and the results were visualized. Finally, a drug sensitivity analysis was carried out to compare differences in sensitivity to various drugs between the high-risk and low-risk groups. Potential drugs were screened using the “pRRophetic” package ([84]http://genemed.uchicago.edu/~pgeeleher/pRRophetic/). All drugs included in the drug sensitivity analysis of this study are listed in Supplementary 4. COAD sample clustering and inter-cluster analysis We performed consensus clustering analysis on COAD samples, dividing them into clusters ranging from 2 to 9 based on the expression levels of prognostic NRGlncRNAs. We then used a consensus cumulative distribution function plot to visually assess and determine the optimal number of clusters for sample clustering. Survival analysis was carried out between different clusters using R, and survival curves were plotted. The “ggalluvial” package was used to draw an alluvial plot to intuitively see the correspondence between different clusters and high-risk and low-risk groups (10.21105/joss.02017). Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were performed on different clusters and high-risk and low-risk groups to reduce dimensionality and determine differences and similarities between clusters and the risk groups. The “estimate” package was used to score stromal cells, immune cells, and ESTIMATE (the sum of the first two) of TME in different cluster samples, to compare TME differences between clusters. Immune cell infiltration analysis was used to explore the differences in immune cell expression between different clusters and a heat map was plotted. Immune checkpoint analysis involved differential analysis of immune checkpoint genes in each cluster and was represented by box plots. Finally, drug sensitivity analysis was performed among clusters for differential comparison, and the results were visualized as box plots. Colon cancer cell line culture Three cell lines were purchased from Procell: human intestinal epithelial cells (NCM460) and two colorectal cancer cell lines (HCT116 and SW620). All cells were cultured in a medium with 10% fetal bovine serum (Gibco BRL, Gaithersburg, MD, USA) at 37 °C, 95% humidity, and 5% CO[2]. Necroptosis induction Necroptosis was induced using the Necroptosis Inducer Kit with TSZ (Beyotime, C1058S), following the manufacturer’s instructions. To collect the conditioned medium, cells were washed twice with phosphate-buffered saline, treated with the reagents for 4 h, then replaced with fresh culture medium and incubated at 37 °C for 12 h. The conditioned medium was then collected, centrifuged at 1500 rpm for 5 min to collect the supernatant, and stored at 4 °C. Quantitative reverse transcription PCR Total RNA was extracted and purified using Trizol reagent (Beyotime, R0016), then reverse transcribed into cDNA. Quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR) was performed using the BeyoFast™ SYBR Green One-Step qRT-PCR Kit (Beyotime, D7268S) to measure relative lncRNA expression using the 2^−ΔΔCt method. Primer sequences used in this study are listed in Supplementary 5. Statistical analysis We employed Strawberry Perl 5.32.1.1([85]https://strawberryperl.com/) for the extraction and annotation of TCGA and GEO datasets, while all other statistical analyses were performed using R V4.2.1([86]https://www.r-project.org/). For two groups, we used a t-test or the Wilcoxon rank-sum test. For data across three or more groups, we used a one-way analysis of variance and the Kruskal-Wallis rank-sum test. Correlation analysis was carried out using the Spearman rank correlation test. We set the P-value at < 0.05 or a false discovery rate (FDR) < 0.05 after correction by the Benjamini-Hochberg method to be statistically significant. Cluster analysis was performed using the k-means algorithm with Euclidean distance and a maximum of nine clusters. In the MR analysis, the random-effects inverse variance weighted (IVW) method was used to assess the causal relationship between prognostic NRGs and COAD; Cochran’s Q statistic was applied for IVW analysis and MR-Egger analysis to detect heterogeneity in the effects of prognostic NRGs on COAD outcomes. A P-value > 0.05 indicated no heterogeneity; MR-Egger regression was used to identify potential horizontal pleiotropy. A P-value > 0.05 indicated no horizontal pleiotropy^[87]36. Results NRGlncRNA analysis By extracting the expression levels of NRGs and combining them with the expression levels of lncRNAs, we filtered according to certain selection criteria and obtained the co-expression relationships between the two. Approximately 45% of NRGs showed a significant co-expression relationship with lncRNAs, and almost all co-expression relationships were positive. The primary results are presented in Fig. [88]2A. Differential analysis was conducted on the NRGlncRNA expression levels in the normal and tumor sample groups, followed by filtering based on certain selection criteria. The expression levels of the top 100 most significantly up-regulated and down-regulated NRGlncRNAs are represented in a heatmap, while the overall results are displayed in a volcano plot, as shown in Fig. [89]2B and C. We can see that the number of up-regulated NRGlncRNAs exceeds that of the down-regulated ones. Fig. 2. [90]Fig. 2 [91]Open in a new tab (A) Co-expression network of NRGs and lncRNAs; (B) Heatmap of expression of NRGlncRNAs; (C) Volcano plot of expression of NRGlncRNAs; (D) Forest plot of hazard ratio of prognostic NRGlncRNAs; (E) Heatmap of expression of prognostic NRGlncRNAs; (F) Alluvial plot of prognostic NRGs and NRGlncRNAs; (G) Heatmap of expression of prognostic NRGs. (*** : P < 0.001, ** : P < 0.01, and * : P < 0.05). Construction of the NRG-lncRNA prognostic model Survival analysis was conducted on clinical information, including survival time, survival status, and NRGlncRNA expression levels, yielding significant NRGlncRNAs and their relationships with the hazard ratio, as shown in Fig. [92]2D. From the figure, we can observe 22 survival-related NRGlncRNAs, of which 21 are high-risk NRGlncRNAs and one is low-risk NRGlncRNAs. Next, the samples were randomly and evenly divided into a train group and a test group. Cox analysis, Lasso regression analysis, and cross-validation were then conducted. We screened out NRGlncRNAs with significant p-values and derived the model formula: graphic file with name M1.gif Table [93]1 presents the information on NRGlncRNAs and coefficients. According to the differential expression analysis of the 22 prognostic NRGlncRNAs in the tumor and normal groups, as shown in Fig. [94]2E, we can observe that apart from LINC0857, [95]AL161785.1, and LINC02381, which have higher expression in the normal group, the remaining NRGlncRNAs are highly expressed in the tumor group. By constructing an alluvial plot (Fig. [96]2F), we identified 11 NRGs co-expressed with the 22 prognostic NRGlncRNAs, all of which have positive regulatory relationships. The differential analysis we conducted on these prognostic NRGs showed that nine of them have statistical differences (AXL, BACH2, CFLAR, CYLD, and IPMK have lower expression in the tumor group, while MAP3K7, ATRX, BRAF, and OTULIN have higher expression), as depicted in Fig. [97]2G. Enrichment analysis of these genes identified their involvement in biological functions such as nucleotide-binding oligomerization domain containing signaling pathway, response to an inorganic substance, regulation of MAPK cascade, and regulation of lymphocyte differentiation. The pathways involved include the nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) signaling pathway. Table 1. The information of NRGlncRNAs and coefficients of the formula of the prognosis model. i NRGlncRNA ID Coefficients 1 AC245100.5 0.576292978 2 [98]AP001619.1 1.408412851 3 LINC01614 0.353244999 4 [99]AC010463.3 0.870212996 5 LINC00513 -1.03003968 6 [100]AL162595.1 0.657806078 7 ITGB1-DT 1.441390546 8 LINC01857 0.890007412 [101]Open in a new tab Validation of prognostic NRGs and NRGlncRNAs using the GEO datasets Differential expression analysis of prognostic NRGs in colon cancer samples and normal intestinal samples from the [102]GSE44076 dataset revealed that, of the 11 prognostic NRGs identified in the TCGA dataset, ten are contained in this dataset, of which nine NRGs showed statistically significant differences, as shown in Fig. [103]3A and B. Specifically, ATRX, BRAF, and MAP3K7 showed higher expression in the tumor group, while AXL, BACH2, CFLAR, CYLD, MAPK8, and TSC1 showed higher expression in the normal group. Furthermore, differential expression analysis of prognostic NRGlncRNAs between COAD samples and adjacent normal tissue samples from the [104]GSE70880 dataset (as shown in Fig. [105]3C and D) indicated that, of the 22 prognostic NRGlncRNAs identified in the TCGA dataset, seven are included in this dataset. Among them, three prognostic NRGlncRNAs showed statistically significant differences, with LINC01655 and LINC01614 showing higher expression in the tumor group and [106]AL161785.1 showing higher expression in the normal group. Fig. 3. [107]Fig. 3 [108]Open in a new tab (A) Boxplot of GEO validation in prognostic NRG expression difference; (B) Heat map of GEO validation in prognostic NRG expression difference; (C) Boxplot of GEO validation in prognostic NRGlncRNA expression difference; (D) Heat map of GEO validation in prognostic NRGlncRNA expression difference. (*** : P < 0.001, ** : P < 0.01, and * : P < 0.05). Results of MR analysis SNP characteristics of the eleven prognostic NRGs can be found in (Supplementary 6—Information of Genetic Instrumental Variants of Prognostic NRGs), and none of the SNPs are weak instrumental variables. In the adjusted results after Bonferroni correction, p < 0.0045 suggests strong evidence, and p within the range of 0.0045 < p < 0.05 suggests preliminary evidence. The IVW method indicated that out of the 11 prognostic NRGs, seven showed positive causal relationships with COAD. Strong evidence supports positive causal relationships for AXL, BACH2, CFLAR, CYLD, and TSC1 with COAD, while preliminary evidence suggests positive causal relationships between IPMK, MAP3K7, and COAD(Supplementary 7—MR Analysis between Prognostic NRGs and COAD). ATRX has zero SNP selected under the instrumental variable conditions, thus an MR analysis of ATRX and COAD is not available. Among the seven prognostic NRGs, AXL, CYLD, and IPMK negatively regulate and inhibit the onset of COAD, while BACH2, CFLAR, MAP3K7, and TSC1 positively regulate and promote the onset of COAD (Supplementary 7—MR Analysis between Prognostic NRGs and COAD). Through heterogeneity analysis, the results indicate the reliability of IVW results and no heterogeneity in MR analysis results of prognostic NRGs and COAD (P > 0.05) (Supplementary 8—Heterogeneity and Horizontal Pleiotropy Results of MR analysis between Prognostic NRGs and COAD). MR-Egger regression results indicated that the MR analysis results for IPMK exhibit horizontal pleiotropy (P = 0.03), thus rendering the MR analysis results for IPMK unreliable (Supplementary 8—Heterogeneity and Horizontal Pleiotropy Results of MR analysis between Prognostic NRGs and COAD). (MR analysis results in Fig. [109]4A and G). Fig. 4. [110]Fig. 4 [111]Open in a new tab MR analysis results of prognostic NRG with COAD. (A) Scatter plot of AXL and COAD using different MR methods. The slopes of the lines represent the causal effects estimated by each method; (B) Scatter plot of BACH2 and COAD using different MR methods; (C) Scatter plot of CFLAR and COAD using different MR methods; (D) Scatter plot of CYLD and COAD using different MR methods; (E) Scatter plot of MAP3K7 and COAD using different MR methods; (F) Scatter plot of TSC1 and COAD using different MR methods; (G) Forest plot of prognostic NGR and COAD using IVW. Prognostic analysis and construction of a nomogram By conducting survival analysis on high-risk and low-risk patients in the train group, test group, and the overall group, statistical differences in the survival curves of high-risk and low-risk patients can be observed across all three groups, as shown in Fig. [112]5A and C. According to the risk curve based on patients’ risk scores, the number of deaths in the train, test, and overall groups increases with the rise in risk score, consistent with our prediction. Simultaneously, we can ascertain that, among the eight prognostic risk model NRGlncRNAs, LINC00513 is lowly expressed while the remaining seven NRGlncRNAs are highly expressed, as depicted in Fig. [113]5D and F. From the integrated results of univariate and multivariate independent prognostic analyses in Fig. [114]5G and H, it can be seen that clinical characteristics such as age, gender, stage, and TNM staging system cannot serve as independent prognostic factors. However, the risk score from our constructed model can act as a prognostic factor independent of its clinical features. From the ROC results in Fig. [115]5I and J, our constructed model exhibits high accuracy and is more accurate than other clinical features. Subsequently, we validated the model across different clinical groups and found that it is applicable to early-stage and late-stage patients, both males and females, as well as patients both above and below 65 years of age, as shown in Fig. [116]5K and P. The generated nomogram (Fig. [117]5Q) predicts the patient’s survival time and prognosis by calculating the total score of each clinical feature. Meanwhile, as can be seen from Fig. [118]5R, the nomogram exhibits high accuracy. Fig. 5. [119]Fig. 5 [120]Open in a new tab (A-C) Survival curves of train group, test group and overall group; (D-F) Heatmap of expression differences of prognostic model NRGlncRNAs in Train group, Test group and Total group; (G-H) Univariate and multivariate independent prognostic analysis; (I) The ROC curves of the risk score and clinical features; (J) ROC curve of 1, 3, 5-year survival rate based on prognostic model NRGlncRNAs; (K-P) Survival curves for different subgroups (ages, genders, and stages); (Q, R) Nomogram of the risk score and clinical features and its calibration plot. Results of GSEA, immune cell correlation analysis, ssGSEA, tumor microenvironment analysis, checkpoint analysis, and drug sensitivity screening Significant pathways enriched in the high-risk group, as screened by GSEA, include the phosphatidylinositol signaling system, transforming growth factor beta (TGF-β) signaling pathway, pathways in cancer, Janus kinase (JAK) - signal transducer and activator of transcription (STAT) signaling pathway, and FC epsilon RI signaling pathway. The low-risk group showed significant enrichment in pathways such as peroxisome, glutathione metabolism, pyruvate metabolism, oxidative phosphorylation, and pentose phosphate pathway, as shown in Fig. [121]6A. Immune cell correlation analysis using seven methods revealed that 12 types of immune cells are significantly negatively correlated with risk scores, while 27 types of immune cells show a significant positive correlation. This indicates that the majority of immune cells are positively correlated with risk scores, as depicted in Fig. [122]6B. Significantly differentially expressed immune functions derived from ssGSEA are all highly expressed in the high-risk group. At the same time, significantly differentially expressed immune cells are all highly expressed in the high-risk group, as shown in Fig. [123]6C and D. From Fig. [124]6E and F, it can be seen that TME analysis reveals statistical differences in stromal cell scores, immune cell scores, and estimate scores between high and low-risk score groups. Immune checkpoint analysis reveals that 36 immune checkpoint genes show statistically significant differences between high and low-risk groups, all of which are highly expressed in the high-risk group, as shown in Fig. [125]6G. Drug sensitivity analysis based on our constructed risk model revealed statistically significant differences in the sensitivity to 54 drugs between high-risk and low-risk groups (Supplementary 9—Drug Sensitivity Analysis). Among these, 47 drugs had significantly lower IC50 values in the high-risk group, suggesting higher sensitivity. Only 7 drugs showed higher sensitivity in the low-risk group, suggesting that most drugs could have potential therapeutic effects on the high-risk population. Fig. 6. [126]Fig. 6 [127]Open in a new tab (A) Curves of GSEA enriched KEGG pathways between high and low risk groups; (B) Bubble plot of immune cell correlation analysis based on seven methods; (C, D) Boxplot of significantly differentially expressed ssGSEA immune functions and immune cells between high and low risk groups; (E, F) Boxplot of the difference in TME score of stroma and immune cells between high and low risk groups; (G) Boxplot of differential expression of immune checkpoint genes between high and risk groups. (*** : P < 0.001, ** : P < 0.01, and * : P < 0.05). Results of COAD sample clustering and inter-cluster analysis We divided the tumor samples into three clusters, namely C1, C2, and C3, through cluster analysis, as shown in Fig. [128]7A. Survival analysis of these three clusters revealed statistical differences between their survival curves of the three clusters, as shown in Fig. [129]7B. The subsequent construction of an alluvial plot (Fig. [130]7C) clarified the relationships between C1, C2, and C3 and high and low-risk scores. PCA and t-SNE analysis (Fig. [131]7D and G) revealed that C1 and C2 exhibited certain similarities, and their main differences were with C3. Specifically, C3 overlapped mainly with the high-risk group, C2 overlapped mainly with the low-risk group, and C1 had a small gap with both. TME analysis of the clustered samples showed that there were statistical differences between C1, C2, and C3 in terms of stromal cell scores. There were also statistical differences between C1 and C3, and C2 and C3 in immune cell scores, as well as estimate scores, as shown in Fig. [132]7H and J. Immune cell correlation analysis showed that C3 mainly exhibited upregulation of immune cells, whereas C1 and C2 mainly exhibited downregulation, as shown in Fig. [133]7K. Immune checkpoint analysis revealed that the expression of 40 immune checkpoint-related genes showed statistical differences. Among these, the expression of 32 immune checkpoint genes was upregulated in C3, the expression of 7 genes was upregulated in C1, and one immune checkpoint gene (TNFSF15) was upregulated in C2, as shown in Fig. [134]7L. Upon conducting a drug sensitivity analysis among clusters, we found that 74 drugs showed sensitivity differences among the three clusters (Supplementary 10—Drug Sensitivity Analysis among Clusters). Of these, 70 drugs had lower IC50 in C3, and 43 drugs corresponded to those identified in the sensitivity analysis between high and low-risk groups. Fig. 7. [135]Fig. 7 [136]Open in a new tab (A) Cluster analysis results of COAD samples; (B) Survival curves among different clusters; (C) Alluvial plot of association between clusters and risk groups; (D-G) PCA scatter plot of high-risk and low-risk groups and clusters; (H-J) Boxplot of the difference in TME score of stroma, immune cells, and their sum among clusters; (K) differential heat map of immune cell expression among different clusters; (L) Boxplot of differential expression of immune checkpoint genes among clusters. (*** : P < 0.001, ** : P < 0.01, and * : P < 0.05). Validation results of prognostic model NRGlncRNAs in cell experiments To validate the function of the eight prognostic model NRGlncRNAs in COAD, we examined the expression levels of these prognostic NRGlncRNAs in the normal intestinal epithelial cell line NCM460 and the colorectal cancer cell lines HCT116, and SW620. The results are detailed in Supplementary 11. The Wilcoxon signed-rank test was used to analyze the significance of the expression differences of NRGlncRNAs between the two colorectal cancer cell lines and NCM460. The results showed that AC245100.5, [137]AP001619.1, LINC01614, [138]AC010463.3, [139]AL162595.1, ITGB1-DT, and LINC01857 were highly expressed in the colorectal cancer cell lines, while LINC00513 was lowly expressed in the colorectal cancer cell lines (P < 0.05, as shown in Fig. [140]8A and B). Fig. 8. Fig. 8 [141]Open in a new tab (A) Differential expression analysis results of prognostic NRGlncRNAs between the tumor group (HCT116) and the normal group (NCM460); (B) Differential expression analysis results of prognostic NRGlncRNAs between the tumor group (SW620) and the normal group (NCM460). (* : P < 0.05) Discussion Necroptosis, as a form of programmed cell death distinct from apoptosis, is closely linked with intestinal diseases, being involved in the occurrence, development, and prognosis of ulcerative colitis (UC), Crohn’s disease and other inflammatory bowel diseases (IBD)^[142]37,[143]38, various intestinal injuries^[144]39–[145]41, colonic adenomas, and colorectal cancer (CRC)^[146]42,[147]43. Necroptosis has a dual role in cancer. It can cause cancer cell death and thus benefit cancer treatment, but it can also cause inflammation leading to tumor progression through immune pathways^[148]44. As the most common intestinal malignancy^[149]45, COAD is linked to necroptosis. Current research suggests that necroptosis can trigger necroptosis-dependent immune escape, thus promoting the progression of colon cancer^[150]46. However, opposing studies argue that necroptosis can suppress inflammatory responses, thus alleviating colitis and reducing the risk of colon cancer^[151]47. Additionally, some research suggests that necroptosis has no significant role in colon cancer^[152]48. Given these conflicting findings, the role of necroptosis in COAD remains unclear. Our study aims to further elucidate the regulatory mechanisms of necroptosis in COAD and explore potential therapeutic strategies. Firstly, we identified 11 prognostic NRGs from the TCGA-COAD dataset: ATRX, AXL, BACH2, BRAF, CFLAR, CYLD, IPMK, MAP3K7, MAPK8, OTULIN, and TSC1. We selected the COAD dataset from the Chinese population in the GEO database for validation. Nine of the prognostic NRGs showed statistically significant differences. The consistency in different population datasets supports the reliability of these study results, indicating that these prognostic NRGs are strongly correlated with COAD. According to the MR analysis results, there are six prognostic NRGs associated with the occurrence and development of COAD, three of which are consistent with the study results of both the TCGA and GEO datasets and satisfy the causal relationship with COAD. Specifically, the expression levels of AXL and CYLD are causally related to a reduced risk of COAD, while the expression level of MAP3K7 is causally related to an increased risk of COAD. CYLD is a tumor suppressor factor that inhibits the proliferation of colon cancer cells by negatively regulating the NF-κB and JNK signaling pathways through interaction with signaling molecules such as NEMO and TRAF2^[153]49–[154]52. The survival of colon cancer cells depends on the activity of MAP3K7, and by inhibiting MAP3K7, the overactivated Wnt signaling pathway can be suppressed, inducing colon cancer cells towards apoptosis; the combined use of MAP3K7 inhibitors and commonly used chemotherapy drugs can effectively reduce tumor drug resistance, thereby inhibiting tumor growth and survival^[155]53,[156]54. AXL is a protein expressed in many tumor cells, and its increased expression is closely related to the enhanced invasion and metastasis of colon cancer cells. Therefore, it is often regarded as a biomarker indicating a poor prognosis of colon cancer^[157]55,[158]56. Current research indicates that AXL is also expressed in macrophages and dendritic cells, playing a key role in suppressing inflammatory responses and reducing the risk of colon cancer^[159]57. Although current anti-cancer therapies targeting AXL are considered a potential strategy, more research is needed to understand its specific role in COAD due to its complex dual role in cancer. BACH2 acts as an intrinsic negative regulator in NK cells, restricting their maturation and function, and thereby inhibiting the metastasis of cancer cells^[160]58. CFLAR is a negative regulator of necroptosis, which can prevent apoptosis signal transduction and regulate other cell death pathways. Due to its aberrant expression in various cancers, CFLAR holds promise as a potential target for cancer therapy^[161]59. IPMK exerts anti-tumor activity by activating the protein kinase AMP-activated catalytic subunit alpha 1 (AMPK) / unc-51 like autophagy activating kinase 1 (ULK1) signaling pathway and inducing autophagy in tumor cells^[162]60. ATRX is one of the most commonly mutated tumor suppressor genes in human cancers. When mutated, it affects cellular maintenance of telomeres and chromatin regulation, thereby influencing cell growth and differentiation, and increasing the risk of cancer^[163]61,[164]62. BRAF mutations lead to sustained MAPK signaling, allowing cancer cells to grow and divide uncontrollably, providing favorable conditions for cancer development. Therefore, BRAF is a crucial target in the treatment of colorectal cancer^[165]63. Phosphorylation of OTULIN promotes genotoxicity-induced Wnt/β-catenin activation, which plays a critical role in cancer progression, metastasis, and drug resistance^[166]64. These prognostic NRGs are mainly involved in regulating ubiquitination status, promoting oxidative stress and mitochondrial dysfunction, and increasing tumor cell starvation induction, thereby modulating necroptosis. Currently, the recruitment and activation of RIPK1, RIPK3, and MLKL are considered critical processes in necroptosis and play a pivotal role in tumor progression^[167]65. Overexpression of AXL and BRAF leads to the loss of RIPK3 expression, enabling cancer cells to escape necroptosis^[168]66. CYLD works synergistically with the necroptotic pathway involving RIPK1, promoting cancer cell apoptosis and exerting anti-tumor effects^[169]67. MAP3K7 inhibits the induction of cell death involving RIPK1, influencing the response of cancer cells to necroptosis and exacerbating cancer progression^[170]68. OTULIN deubiquitinates LUBAC (linear ubiquitin chain assembly complex), suppressing TNF receptor superfamily member 1 A (TNFR1) and RIPK1-mediated necroptosis, thereby limiting cell death and inflammatory responses^[171]69. Due to limited research on lncRNAs, no studies have reported on the relationship between the identified prognostic NRGlncRNAs and the classical necroptosis pathway of RIPK3/RIPK1/MLKL. Further in-depth research is required to explore this association. Enrichment analysis of prognosis NRGs indicates their predominant involvement in the NF-κB signaling pathway. Experimental validation affirms that the phosphorylation activating this pathway plays a mediating role in promoting necroptosis^[172]70. Concurrently, this pathway is associated with the release of inflammatory factors, regulating inflammation via the classical necroptosis pathway involving RIPK3/RIPK1/MLKL^[173]71. Current research confirms that the expression of NF-κB in COAD can mediate tumor progression and metastasis^[174]72. Based on our comprehensive results, it can be inferred that regulators like RIPK3/RIPK1/MLKL may positively modulate the NF-κB pathway, thereby promoting COAD. We identified ten immune cells with statistically significant differences through immune cell analysis. Among them, activated dendritic cells (aDCs) recognize antigens on the surface of tumor cells and stimulate the immune system to attack tumors by presenting these antigens. However, considering the inhibitory mechanisms within tumors, aDCs are often insufficient to induce a potent immune response^[175]73,[176]74. Immune B cells can impact anti-tumor immunity by expressing leucine-tRNA synthetase-2, leading to immune evasion in colorectal cancer cells^[177]75. Macrophages can regulate the infiltration and epithelial-mesenchymal transition program of colorectal cancer cells by secreting interleukin-6, enhancing the migration and invasion capabilities of colorectal cancer cells^[178]76. Neutrophils, through the secretion of anterior gradient 2, protein disulfide isomerase family member (AGR2) acting on its receptor CD98hc-xCT, can promote the metastasis of colorectal cancer^[179]77. NK cells’ immune infiltration can directly attack and destroy colorectal cancer cells, but colorectal cancer cells can counteract immune infiltration by upregulating major histocompatibility complex, class I, E (HLA-E) (ligand of the inhibitory receptor killer cells lectin-like receptor C1 (NKG2A) expressed by CD8 and NK cells), inhibiting the immune activity of CD8 and NK cells^[180]78. Helper T cells can induce cytokine response-related colonic inflammation, enhancing COAD formation^[181]79,[182]80. Follicular helper T cells exhibit their activation role in promoting anti-intestinal tumor immunity^[183]81. Tumor-infiltrating lymphocytes have been confirmed to promote the development of immune heterogeneous tumor microenvironments in colorectal cancer, thereby driving tumor progression^[184]82. The positive correlation between Th1 cells and the survival rate of cancer patients suggests that inhibiting Th1 cells may be a crucial carcinogenic mechanism^[185]83. Modulating T cell recruitment can reduce anti-tumor CD8 T cell immunity, thus promoting colorectal cancer^[186]84. The specific expression characteristics of lncRNA play a key role in cancer progression and survival prediction^[187]85. We employed Lasso and Cox regression analysis to identify eight prognostic NRGlncRNAs and validate them using the GEO dataset and cell experiments. One of the NRG-lncRNAs in the GEO dataset, LINC01614, was consistent with the results found in the TCGA dataset. The cell experiment results also aligned with those from the TCGA dataset. Of the eight prognostic NRGlncRNAs, seven were highly expressed in tumor tissues, namely AC245100.5, [188]AP001619.1, LINC01614, [189]AC010463.3, [190]AL162595.1, ITGB1-DT, LINC01857, while LINC00513 was lowly expressed in tumor tissues. Research on the role of these eight prognostic NRGlncRNAs in COAD is still limited. LINC01614 is an oncogenic lncRNA that directly interacts with ANXA2 and p65, promoting the activation of NF-κB, resulting in the upregulation of glutamine transporter proteins SLC38A2 and SLC7A5, thereby facilitating the malignant progression of cancer^[191]86. High expression of ITGB1-DT can enhance the proliferation, migration, and invasion of lung adenocarcinoma cells. It is significantly associated with a later clinical stage of lung adenocarcinoma, a lower overall survival rate, and a lower disease-free survival rate^[192]87. LINC01857 is also an oncogenic lncRNA that promotes the growth, migration, and invasion of cancer cells by modulating the miR-1281/TRIM65 axis^[193]88. LINC00513 is a regulator of the interferon signaling pathway, influencing the tumor microenvironment by modulating inflammatory responses^[194]88. Although the biological functions of AC245100.5, [195]AP001619.1, [196]AC010463.3, and [197]AL162595.1 have not been confirmed in relevant research, their potential mechanisms linking these NRGlncRNAs to the pathogenesis of COAD are of significant research value. We performed GSEA on the eight NRGlncRNAs used to construct the prognostic model and found that they are mainly involved in the following pathways: phosphatidylinositol signaling system, TGF-β signaling pathway, cancer signaling pathway, JAK-STAT signaling pathway, FC epsilon RI signaling pathway, peroxisome, glutathione metabolism, pyruvate metabolism, oxidative phosphorylation, and pentose phosphate pathway. Current research indicates that the phosphatidylinositol-related Phosphoinositide 3-kinase (PI3K) / threonine kinase 1 (AKT) signaling pathway and the TGF-β signaling pathway have carcinogenic potential and can promote the progression of CRC^[198]90–[199]92. The activation of the JAK-STAT pathway can lead to the proliferation and migration of colon cancer cells^[200]93. The exact function of FC εRI in COAD has not been confirmed, but its expression is found in patients with COAD and gastrointestinal inflammation^[201]94. Peroxisome proliferator-activated receptor γ (PPARγ) can downregulate cell proliferation in colon cancer cells^[202]95. Elevated levels of serum glutathione and enhanced activity of glutathione reductase and glutathione-S-transferase can inhibit oxidative stress, reduce inflammation, reduce cell proliferation, and increase cell apoptosis to prevent colon cancer^[203]96. The upregulation of pyruvate kinase M2 (PKM2) can have a therapeutic effect on CRC^[204]97. Oxidative phosphorylation can weaken the proliferation of CRC cells and the growth of CRC tumors^[205]98. Enzymes in the pentose phosphate pathway can increase DNA synthesis and proliferation in CRC cells, thus promoting tumor progression^[206]99. Furthermore, immune cell correlation analysis, ssGSEA, tumor microenvironment analysis, and immune checkpoint analysis all suggest the same result: the expression of immune cells, immune functions, tumor microenvironment, and immune checkpoint genes are positively correlated with risk score. This high level of consistency suggests that NRG-related immune inflammation in COAD primarily promotes tumor progression. Finally, through drug sensitivity analysis, we found that the high-risk group in our risk model mostly has lower IC50 values for drugs, suggesting a higher sensitivity to treatment. We identified eight NRGlncRNAs for constructing the prognostic model that can modify prognostic NRGs and mediate signaling pathways such as NF-κB, JAK-STAT, and PI3K/AKT. These NRGlncRNAs also regulate various immune cells, including T cells, B cells, and macrophages. Ultimately, they predominantly promote COAD development through the necroptotic process. The biological functions involved encompass immune inflammatory response, oxidative stress, immune escape, telomere regulation, and cytokine response. In addition, this study employed the Lasso + Cox method, which is widely used in constructing prognostic models for gene expression data in biomedical research. By maximizing the elastic net penalized partial likelihood, this method offers better stability and computational speed compared to other competing algorithms^[207]100. Deep learning methods are primarily applied in the biomedical field for medical image analysis, diagnostic model construction, and biomedical text mining^[208]101,[209]102. These methods exhibit high accuracy and are easy to train, but they still require optimization for handling imbalanced datasets and rely on substantial computational power. Nevertheless, we believe that deep learning has broad application prospects in gene expression data analysis. Conclusion This study systematically analyzed the role of prognostic NRGlncRNAs in predicting the survival of COAD patients and their functions in the tumor microenvironment, immune checkpoints, and immune cell infiltration. Additionally, this study explored potential regulatory mechanisms of prognostic NRGlncRNAs and identified potential drugs for COAD treatment. The study identified eight NRGlncRNAs for constructing the prognostic model and nine prognostic NRGs with regulatory relationships to them. The findings were analyzed and validated using GEO and GWAS datasets, as well as cell experiments. These NRGlncRNAs and NRGs can predict the survival rate of COAD patients, providing clinical guidance for the treatment of COAD patients. Although this study employed multiple datasets for multidimensional bioinformatic analysis to comprehensively understand the role of necroptosis in COAD, identify prognostic NRGlncRNAs, construct an accurate prognostic risk model, and conduct multifaceted validation of the results, there are still some limitations. This study conducted validation at the cellular level but lacks validation in clinical samples, and further evaluation is required in clinical settings to assess its stability and practical utility. Additionally, although multiple datasets were selected for analysis, the chosen datasets may present sample bias, particularly in relation to the clinical background of patients, racial diversity, and other factors that may influence the progression of COAD. Therefore, future research should consider a broader and more diverse patient population to enhance the generalizability and clinical applicability of the findings. Electronic supplementary material Below is the link to the electronic supplementary material. [210]Supplementary Material 1^ (1.9MB, pdf) Acknowledgements