We downloaded the RNA-sequencing and clinical data of 328 tissue samples, including 40 adjacent normal tissues and 288 CRC tissues, from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/ccg/research/genome-sequencing/tcga). The robust multi-array averaging (RMA) algorithm in the ‘affy’ R package (version 1.46.1) [11] was used to process these data. In addition, batch effects were removed using the ‘sva’ R package (version 3.42.0) [12]. Differentially expressed genes (DEGs) between CRC and normal tissues were identified using the ‘limma’ R package (version 3.42.2) [13], with the screening criteria set as $|$logFC$|$ values of $\ge$1 (FC: fold change) and corrected p-values of 0.05. NRGs were collected from the GeneCards database (https://www.genecards.org/), and a Venn diagram was generated to represent the intersection between DEGs and NRGs. Genes with correlation coefficients of $>$0.2 were used for subsequent analysis.

The ‘ConensusClusterPlus’ R package (version 1.66.0) was used for consensus clustering analysis of NRGs [14]. Based on the expression of NRGs, two molecular subtypes were generated through K-means clustering. Gene set variation analysis (GSVA) was performed on gene sets extracted from MSigDB (C2.Cp.ke.v7.2) [15] to examine differences in biological functions between the two molecular subtypes.

We compared the prognosis and clinical features of patients with CRC between the two molecular subtypes. To compare overall survival (OS) between the two subtypes, Kaplan-Meier curves were generated using the ‘survival’ and ‘survminer’ R packages (version 4.2.3). The CIBERSORT algorithm was used to estimate the immune scores of 22 immune cell types in each CRC sample based on cell-specific gene signatures [16].

The ‘limma’ R package (version 3.42.2) was used to screen for DEGs between the two molecular subtypes [13], with the screening criteria set as $|$logFC$|$ values of $\ge$1 and adjusted p-values of 0.05. Patients with CRC were randomly assigned to training (n = 287) and test (n = 41) sets, which were used to develop a prognostic risk model [17]. Lasso–Cox regression analysis was implemented to reduce the risk of overfitting and to establish a risk model based on prognosis-associated NRGs. The prognostic NRG-Score was calculated based on the risk model as follows [18]: ${\sum }_i^{n}Coefi\times Ai$. In the aforementioned equation, Coefi represents the risk coefficient and Expi represents the expression of the gene. The training and test sets were divided into high- and low-risk groups based on the median NRG-Score. Kaplan-Meier and receiver operating characteristic (ROC) curves were plotted to evaluate the generalisability and predictive accuracy of the model.

The ‘RMS’ package (version 3.6.1) was used to establish a predictive nomogram based on NRG-Scores and clinical characteristics of CRC [19]. Column charts demonstrating NRG-Scores and predicted OS rates were generated, and calibration curves were plotted to compare the predicted and actual 1-, 3- and 5-year OS rates [20].

The infiltration levels of 22 types of immune cells were evaluated based on the expression of the four prognostic NRGs included in the risk model (model genes). A boxplot was generated to compare the expression of immune checkpoint genes between the low- and high-risk groups. In addition, drug sensitivity was compared between the two risk groups. The ‘pRRophetic’ package (version 4.3.2) was used to calculate the half-maximal inhibitory concentration (IC50) of targeted drugs to evaluate their clinical efficacy in the two risk groups [21].

The ‘Seurat’ package (version 5.0.1) was used to generate a Seurat object and remove poor-quality cells. A standard procedure was used to pre-process the data, and the percentages of gene count, cell count and mitochondrial content were subsequently calculated. Cells with 200 or $>$6000 genes and with $>$5% of transcripts derived from mitochondria were excluded. To normalise the data of each cell, UMI counts were scaled using a scale factor of 10,000. After the logarithmic transformation of the data, the ScaleData function was used to analyse corrected normalised data metrics. Principal component analysis (PCA) was performed using the top 10 variable genes, and the top 5 principal components were visualised and clustered via UMAP (or TSNE). Cells were clustered using the FindClusters function implemented in the Seurat R package.

The human CRC cell lines HCT116, HT29 and HCT15 and the normal human colorectal epithelial cell line NCM460 were purchased from the Cell Center of the Institute of Basic Medical Sciences, Chinese Academy of Medical Sciences. All cells were cultured in DMEM (Gibco, New York, NY, USA) supplemented with 10% foetal bovine serum (FBS) (VivaCell, Shanghai, China) and 1% penicillin-streptomycin (Gibco, New York, NY, USA) in a humidified environment with 5% CO${}_2$ at 37 ℃. The cells were transfected with siRNAs using Lipo8000TM transfection reagent (Biyuntian, Shanghai, China) according to the instruction manual. All cell lines were validated by short tandem repeat (STR) profiling and tested negative for mycoplasma. The target sequences of H2BC18 siRNA were shown in Table 1.

For quantitative reverse transcription polymerase chain reaction (qRT-PCR), total RNA was extracted from CRC cells using Trizol reagent (Takara, Kyoto, Japan) and was reverse transcribed using Prime Script RT Master Mix (Takara, Kyoto, Japan). Quantitative PCR (qPCR) was performed using total 2 µL of mRNA, specific primers and SYBR Premix Ex Taq II (Takara, Kyoto, Japan). The mRNA expression of target genes was normalised to that of GAPDH (internal control) and quantified using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{Ct}}$ method [22]. Primers were synthesised by Sangon Biotech (Sangon, Shanghai, China). The gene-specific primers sequences were shown in Table 2.

The cells were digested utilising 0.25% trypsin and centrifuged at a speed of 1200 rpm/min for 5 min. The supernatant was discarded, and the cell precipitate was resuspended in 1 mL of a complete medium. Subsequently, the cells were counted using a Bovine Bow counting plate (QiuJing, Shanghai, China) with a filled cell. The cell suspension was diluted to a concentration of 1 $\times$ 10${}^3$ cells/well, and 200 µL of the cell suspension was added to each well of a 96-well plate. The cells were maintained in an incubator, and the original medium was replaced with 100 µL of a serum-free medium containing 10% CCK-8 reagent on days 0, 1 and 2. The CCK-8-treated cells were then incubated at 37 ℃ for two hours. Following this, the optical density (OD) of every well was measured using an enzyme marker (MCE, South Brunswick Township, NJ, USA) at 450 nm.

Migration assay: Cells were cultured in a serum-free medium overnight. The following day, the cells were digested using 0.25% trypsin and centrifuged at 1200 rpm/min for 5 min. The supernatant was discarded, and the cell precipitate was washed thrice with PBS. The cells were resuspended in 1 mL of a basal medium and counted on Oxbow counting plates. The concentration of HCT15 cells was adjusted to 4 $\times$ 10${}^4$ cells/well. A total of 600 µL of a medium containing 10% FBS was added to a 24-well plate serving as the lower chamber, and a Transwell insert serving as the upper chamber was placed into the plate. HCT15 cells were added to the upper chamber (200 µL/well) and cultured routinely for 48 h. After the media in the Transwell chamber and plate were discarded, the migrated cells were washed twice with PBS, fixed with 4% paraformaldehyde for 30 min and stained with crystal violet for 15 min. A moistened cotton swab was used to gently remove the cells from the upper chamber. Thereafter, the cells were photographed using a Nikon eclipse80i confocal fluorescence microscope (Nikon, Tokyo, Japan).

Invasion assay: The basal medium was used to dilute a substrate gel at a ratio of 1:8. A total of 60 µL of the diluted gel was added to the chambers, followed by incubation at 37 °C for 2 h. The subsequent steps were the same as those for the migration assay.

Wound healing assay: CRC cells were cultured in 6-well plates until $\ge$90% confluence was achieved. The surface of the cell monolayer in each well was scratched with a sterile 200-µL pipette tip. Floating cells were washed with PBS, and a serum-free medium was added to each well. Photographs were captured at 0 h, 48 h and 72 h after the scratches were made, and changes in the area of the scratches were quantified using the ImageJ software (Image J2x 2.1.5.0, National Institutes of Health, New York, NY, USA) [23].

Total proteins were extracted from cells and CRC tissues using RIPA buffer (Solarbio, Beijing, China). The extracted proteins were quantified via BCA assay, separated on 8–12% sodium dodecyl sulphate–polyacrylamide gels and transferred to polyvinylidene difluoride membranes. The membranes were initially incubated with anti-H2BC18 (1:1000, Boiss, Beijing, China) and anti-GAPDH (1:5000, Bioworld, Nanjing, China) primary antibodies and subsequently incubated with horseradish peroxidase (HRP)-conjugated anti-rabbit IgG antibody (secondary antibody, 1:5000, Bioworld, Nanjing, China). Finally, protein bands were visualised using an ECL reagent (Bioworld, Nanjing, China).

The expression of H2BC18 in clinical tissues was assessed via immunohistochemical (IHC) analysis using anti-H2BC18 antibody (1:500, Boiss, Beijing, China), Paraffin-embedded tissue blocks were cut into 4-µm-thick sections, dewaxed, rehydrated using xylene I/II/III and anhydrous ethanol I/II/III (absolute ethyl alcohol) for 15 minutes each and washed with PBS two times for 5 minutes each. For antigen retrieval, tissue slides were placed in a sectioning rack and slowly immersed into an antigen repair solution. The solution was boiled on high for 5 min and on low for 20 min. Subsequently, the sections were washed with PBS two times for 10 minutes each, blocked with 3% bovine serum albumin (BSA) for 30 minutes at room temperature and incubated with primary antibodies overnight at 4 °C in a humidified incubator. The following day, the sections were washed with PBS (3 $\times$ 5 minutes) and incubated with HRP-conjugated anti-rabbit IgG antibody for 30 minutes at room temperature. Thereafter, the sections were washed with PBS (3 $\times$ 5 minutes) and stained with diaminobenzidine (DAB). Finally, the sections were observed under a white light microscope, and integrated OD (IOD) was calculated using the ImageJ software [23]. Data was expressed as the average OD (AOD), which was calculated as follows: IOD/area.

The R (version 4.2.0) software (New York, NY, USA) was used for statistical analysis and visualisation of results. Student’s t-test was used to compare normally distributed quantitative data and Wilcoxon test was used to compare non-normally distributed quantitative data between groups. Statistical significance was denoted as follows: *, p 0.05; **, p 0.01; ***, p 0.001.

The gene expression data of 328 tissue samples, including 40 adjacent normal tissues and 288 CRC tissues, were extracted from TCGA database. A total of 1831 DEGs were identified between CRC and adjacent normal tissues (Fig. 1A); of which, 945 genes were upregulated and 886 genes were downregulated (Fig. 1B). These DEGs were intersected with 160 NRGs obtained from the GeneCards database using a Venn diagram, and 18 differentially expressed NRGs associated with CRC were eventually obtained (Fig. 1C). Subsequently, a consensus clustering algorithm was used to classify patients with CRC based on the expression of the 18 NRGs. With the optimal k value of 2, the patients were divided into cluster A and cluster B (Fig. 1D). Clinical characteristics (Fig. 1E) and NRG function were compared between the two clusters using Kaplan-Meier (K-M) analysis, immune infiltration analysis and GSVA. K-M curves showed that patients in cluster B had shorter OS than patients in cluster A (Fig. 1F). GSVA showed that cluster A was enriched in pathways related to inhibition of tumour development, including type I interferon signalling pathway and activation of mitochondrial autophagy, whereas cluster B was significantly enriched in tumour-promoting pathways such as DNA replication and chromatin aggregation (Fig. 1G). The infiltration levels of 22 types of immune cells differed significantly between the two clusters, with the infiltration levels of activated CD4 T cells and type 2 T helper cells (Th2) being higher in cluster A than in cluster B (Fig. 1H). These results suggested that the two clusters had different TME-associated properties and NRGs prolonged the OS of patients in cluster A by promoting anti-tumour immunity.

Fig. 1.

Screening and categorisation of overall survival (OS)-related necroptosis-related genes (NRGs) in colorectal cancer (CRC). (A) Heatmap showing the top 60 differentially expressed genes (DEGs) (n = 1813) between CRC and adjacent normal tissue samples. (B) Volcano plot showing up-regulated (n = 945) and down-regulated genes (n = 886). (C) Venn diagram demonstrating the intersection between DEGs and necroptosis-related genes (n = 160). (D) Consensus matrix heatmap defining two clusters and associated regions. (E) Relationship of clusters A and B with clinicopathologic features and NRGs. (F) Kaplan-Meier curves (log-rank test, p 0.001) of OS for clusters A and B. (G) GSVA of NRGs in clusters A and B. (H) Differences in the infiltration levels of 22 types of immune cells between the two clusters. Asterisks indicate p-values (*, p 0.05; **, p 0.01). GSVA, Gene Set Variation Analysis.

Although the two molecular subtypes of necroptosis were found to have different prognostic and immune infiltration patterns, these findings are only applicable to patient populations and are not accurate for assessing the impact of NRGs on the prognosis of CRC. Therefore, we established a prognostic model based on the differentially expressed NRGs between clusters A and B for diagnosing CRC and guiding treatment. Briefly, we randomly divided patients with CRC into training (n = 287) and test (n = 41) sets and identified four NRGs significantly associated with CRC prognosis, namely, CABPB1-IT1, H2BC18, HSPA1L and MIR503HG, via LASSO regression analysis. The NRG-Score was calculated based on these four NRGs as follows: GABPB1-IT1 expression level $\times$ 0.00260 + H2BC18 expression level $\times$ 0.01596 + HSPA1L expression level $\times$ 0.02800 + MIR503HG expression level $\times$ 0.01304 (Fig. 2A,B). Based on the median NRG-Score, patients with CRC were divided into low-risk (n = 134) and high-risk (n = 149) groups. NRGs such as CAMK2A, CASP1, CHMP4A, FTL, H2AC20, H2AC6, H2AC8, PLA2G4C, RIPK1, STAT5A/B and TRADD were differentially expressed between the two groups (Fig. 2C). K-M curves showed that the prediction of prognosis was better in the low-risk group than in the high-risk group (Fig. 2D). Similar results were observed in the test set and the total-sample dataset (Fig. 2E,F). Subsequently, ROC curves were plotted to assess the predictive ability of the model. The results showed that the area under the curve (AUC) values for predicting OS at 1, 3 and 5 years were 0.641, 0.742 and 0.751, respectively, indicating that the model had excellent predictive ability, especially for 3- and 5-year survival (Fig. 2G–I). To assess the clinical applicability of the NRG-Score, we integrated NRG-Scores and clinical characteristics such as age, sex and ISS stage to establish a nomogram to predict OS at 1, 3 and 5 years (Fig. 3A). The calibration curve showed that the predicted OS rates were highly consistent with the actual OS rates, indicating that the nomogram had superior predictive ability (Fig. 3B).

Fig. 2.

Construction and evaluation of a prognostic model. (A) 1000-fold cross-validation for the selection of LASSO regression variables. (B) LASSO regression coefficients of necroptosis-related genes. Each curve corresponds to a gene involved in necroptosis. (C) Expression of 18 OS-associated NRGs in the high- and low-risk groups. (D) Kaplan-Meier (K-M) survival curves for the low- and high-risk groups in the training set. (E) K-M survival curves for the low- and high-risk groups in the test set. (F) K-M survival curves for the low- and high-risk groups in the total-sample dataset. (G) Receiver operating characteristic (ROC) curves for predicting 1-, 3- and 5-year OS in the training set. (H) Receiver operating characteristic (ROC) curves for predicting 1-, 3- and 5-year OS in the test set. (I) Receiver operating characteristic (ROC) curves for predicting 1-, 3- and 5-year OS in the total-sample dataset. (*, p 0.05; **, p 0.01; ***, p 0.001). LASSO, least absolute shrinkage and selection operator.

Fig. 3.

Development and validation of the nomogram. (A) Nomogram integrating NRG-Scores and clinical characteristics such as age, sex, TNM and stage for predicting 1-, 3- and 5-year OS. (B) Calibration curve demonstrating the accuracy of the nomogram in predicting survival at 1, 3 and 5 years. The dashed line represents the performance of an ideal nomogram, whereas the solid green, blue and red lines indicate the performance of the established nomogram. **, p 0.01; ***, p 0.001. TNM, Tumor Node Metastasis.

Studies have shown that necroptosis can promote cell death by enhancing tumour immunogenicity and that activation of the necroptotic factor RIPK1/RIPK3 can lead to upregulation of inflammatory chemokines in the TME, promoting immune cell activation [24, 25, 26]. Therefore, we investigated the association between the four NRGs and immune cell infiltration using the CIBERSORT algorithm. The four NRGs showed a significant correlation with most immune cells (Fig. 4A). In particular, a negative correlation was observed between the expression of NRGs and the infiltration levels of activated NK cells and CD4${}^{+}$ memory T cells (Fig. 4B,C). These results suggested that the NRG-Score was markedly associated with the presence of TILs in the microenvironment of CRC and high-risk patients were susceptible to the immunosuppressive effects caused by the reduced abundance of activated NK and CD4 memory T cells in the TME. Furthermore, the expression of most immune checkpoint genes, including CD274, CD70, HHLA2 and TNFSF9, was significantly higher in the low-risk group, and that of the leukocyte-associated immunoglobulin-producing gene LAIR1 was significantly higher in the high-risk group (Fig. 4D). These results suggest that low-risk patients may respond better to LAIR1-targeted therapies, whereas high-risk patients may respond better to CD274- and CD70-targeted therapies.

Fig. 4.

Immune characterisation of prognostic NRGs in the CRC microenvironment. (A) Correlation between the expression of the four NRGs and the abundance of tumour-infiltrating immune cells. (B) Correlation between the abundance of activated Nature Killer (NK) cells and NRG-Scores. (C) Correlation between the abundance of CD4${}^{+}$ memory T cells and NRG-Scores. (D) Histogram demonstrating differences in the expression of immune checkpoints between the high- and low-risk groups. Yellow and blue columns indicate the low- and high-risk groups, respectively (*, p 0.05; **, p 0.01; ***, p 0.001).

To examine the ability of NRG-Scores to guide the clinical treatment of CRC, we assessed the sensitivity of patients in the high- and low-risk groups to potential targeted drugs. The results showed that patients in the low-risk group responded better to LCK inhibitors (Fig. 5A), third-generation AKT inhibitors (Fig. 5B), JNK inhibitors (Fig. 5C) and the third-generation ABL inhibitor ponatinib (Fig. 5D), whereas patients in the high-risk group responded better to oral PARP inhibitors (Fig. 5E). Altogether, these findings suggest that the abovementioned drugs have potential therapeutic value in the treatment of CRC.

Fig. 5.

Drug sensitivity analysis in the high- and low-risk groups. (A) Differences in the expression of LCK inhibitors between the high- and low-risk groups. (B) Differences in the expression of AKT inhibitors between the high- and low-risk groups. (C) Differences in the expression of JNK inhibitors between the high- and low-risk groups. (D) Differences in the expression of the third-generation ABL inhibitor ponatinib between the high- and low-risk groups. (E) Differences in the expression of PARP inhibitors between the high- and low-risk groups.

To assess the expression of prognosis-associated NRGs in the TME of CRC at the single-cell level, we extracted scRNA-sequencing data from the GSE4158911 and GSE4158912 datasets. After quality control and filtering, cells were classified as mast cells, fibroblasts, endothelial cells, epithelial cells and myeloid cells through dimensionality reduction (Fig. 6A,B). HSPA1L was highly expressed in fibroblasts and endothelial cells, whereas MIR503HG was highly expressed in epithelial cells. Given that CRC is a malignant tumour originating from epithelial cells, upregulated NRGs in epithelial cells may promote the transformation of normal epithelial cells to malignant tumour cells, suggesting that NRGs play an important role in the development of CRC (Fig. 6C).

Fig. 6.

Analysis of the expression of NRGs in cell clusters using scRNA-sequencing data. (A) Cell clusters annotated in the GSE4158911 dataset. (B) Cell clusters annotated in the GSE4158912 dataset. (C) Differential expression of two prognostically relevant NRGs in six cell clusters.

To verify the differential expression and molecular functions of NRGs in CRC, we assessed the mRNA expression of the four prognosis-related NRGs (GABPB1-IT1, H2BC18, HSPA1L and MIR503HG) in three CRC cell lines (HCT15, HCT116 and HT29) and a normal colorectal epithelial cell line (NCM460). The expression of GABPB1-IT1 (Fig. 7A), HSPA1L (Fig. 7B) and MIR503HG (Fig. 7C) was significantly lower in HCT15 and HCT116 cells than in NCM460 cells but was significantly higher in HT29 cells than in NCM460 cells. Only the expression of H2BC18 was significantly higher in all three CRC cell lines than in the normal control cell line (Fig. 7D). Consistently, western blotting showed that the protein expression of H2BC18 was significantly elevated in the three CRC cell lines (Fig. 7E). Clinical tissue samples, including 3 CRC and colorectal inflammation tissues, were collected to verify these results. Immunohistochemical (IHC) analysis showed that the protein expression of H2BC18 was higher in CRC tissues than in colorectal inflammation tissues (Fig. 7F). In addition, IHC data from the Human Protein Atls (HPA) database validated the differential expression of H2BC18 protein between CRC and colorectal inflammation tissues (Fig. 7G). Altogether, these results suggest that H2BC18 plays an important role as a prognosis-associated gene in CRC. Therefore, in subsequent experiments, we examined the biological function of H2BC18 in CRC progression.

Fig. 7.

Expression of prognosis-related NRGs in colorectal cancer (CRC) cell lines. (A) Expression of GABPB1-IT1 in CRC cells versus NCM460 cells. (B) Expression of HSPA1L in CRC cells versus NCM460 cells. (C) Expression of MIR503HG in CRC cells versus NCM460 cells. (D) Expression of H2BC18 in CRC cells versus NCM460 cells (*, p 0.05; **, p 0.01; ****, p 0.0001; ns means the difference between the two groups was not statistically significant. Each experiment was repeated three times). (E) Western blotting was performed to evaluate the protein expression of H2BC18 in three CRC cell lines. (F) Comparison of H2BC18 expression between CRC tissues and colorectal inflammation tissues via immunohistochemical (IHC) analysis. (G) Comparison of H2BC18 expression between CRC tissues and adjacent normal tissues using IHC data from the HPA database. Tumor image available from: https://images.proteinatlas.org/7814/18508_A_8_3.jpg and Normal image available from: https://images.proteinatlas.org/48671/113841_A_7_3.jpg.

To assess the effects of H2BC18 on the malignant progression of CRC, HCT15 cells were transiently transfected with an siRNA targeting H2BC18. The proliferative, invasive and migratory abilities of CRC cells were examined after the successful knockdown of H2BC18 (Fig. 8A). CCK-8 assay showed that knockdown of H2BC18 inhibited the proliferation of HCT15 cells (Fig. 8B). In addition, wound healing and transwell invasion assays showed that knockdown of H2BC18 significantly inhibited the migratory and invasive abilities of HCT15 cells (Fig. 8B–D). Altogether, these results suggest that H2BC18 plays an important role in regulating the malignant progression of CRC.

Fig. 8.

H2BC18 regulates the proliferative, invasive and migratory abilities of CRC cells. (A) Western blotting was performed to detect the efficiency of H2BC18 knockdown in HCT15 cells. (B) CCK-8 assay was performed to examine the proliferative ability (OD450) of CRC cells after H2BC18 knockdown. (C) Wound healing assay was performed to examine the migratory ability of CRC cells after H2BC18 knockdown. (D) Transwell assay was performed to examine the invasive and migratory abilities of CRC cells after H2BC18 knockdown (*, p 0.05; **, p 0.01; ****, p 0.0001).