Introduction
Colorectal cancer (CRC) is a major cause of cancer deaths globally, being the third most diagnosed and second deadliest cancer.1 Its poor prognosis is mainly due to treatment resistance. Despite advances in surgery, chemotherapy, and targeted therapies, survival rates for advanced CRC remain low.2 Immune checkpoint inhibitors (ICIs) like pembrolizumab and nivolumab show promise for metastatic CRC with mismatch-repair deficiency (dMMR) and high microsatellite instability (MSI-H), offering potential long-term remission for these patients.3 However, treating CRC with low immunogenicity, such as MSI-L/pMMR tumors, remains challenging. Current research aims to refine patient selection and develop new combination strategies to boost the effectiveness of immune checkpoint inhibitors.4 Consequently, more studies are urgently needed to identify biomarkers that predict the success of targeted and immune therapies in CRC and to explore their underlying mechanisms.
Neddylation is a post-translational modification where NEDD8 is attached to specific proteins via a three-step enzymatic process involving E1, E2, and E3 enzymes.5 It targets cullin and non-cullin proteins, activating cullin RING ligases and affecting the stability and function of non-cullin proteins.6 This pathway is often overactive in cancers.7–9 Inhibiting neddylation with the NAE inhibitor MLN4924 (pevonedistat) shows strong anticancer effects by promoting cell death10–13 and enhancing cancer cell sensitivity to various treatments.14–18
Beyond its impact on cancer cells, neddylation significantly influences the tumor microenvironment (TME) by modulating the activities of immune cells, such as macrophages, dendritic cells (DCs), and T cells, as well as cancer-associated fibroblasts (CAFs), cancer-associated endothelial cells (CAEs), and various other factors.19,20 The inhibition of neddylation results in the upregulation of PD-L1 expression and the suppression of cancer-associated immunity,21 highlighting neddylation as a promising target for cancer therapy. Moreover, neddylation-related patterns and scores offer potential for distinguishing TME characteristics, predicting patient prognosis, and guiding personalized treatment strategies in oncology.22,23 However, the specific effects of neddylation and neddylation-related genes (NRGs) on CRC remain not well-studied.
Bulk RNA sequencing (RNA-seq) is extensively utilized in translational research to ascertain the average transcript expression within heterogeneous cell populations.24 In contrast, single-cell RNA sequencing (scRNA-seq), a more recent advancement, evaluates gene expression at the individual cell level, thereby elucidating the distribution of expression across various cell subtypes.25 The application of machine learning techniques enhances the analysis of data derived from both methodologies, thereby refining insights into diseases such as cancer and facilitating their clinical application.26–28
In this study, we conducted a comprehensive mapping of neddylation-related genes (NRGs) across 16 distinct cell types in colorectal cancer (CRC) at the single-cell resolution. We performed an analysis of functional variations in neddylation patterns and employed ten machine learning algorithms to develop a gene signature, referred to as NRGS. This signature is designed to predict patient prognosis, characterize the tumor microenvironment, and assess responses to immunotherapy and chemotherapy, utilizing bulk RNA sequencing data.
Methods
Data Acquisition
Single-cell transcriptome datasets GSE166555 and GSE139555 were accessed from the Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org/home/). Bulk RNA-seq data were obtained from the UCSC Xena browser (https://xenabrowser.net/) and the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The TCGA-COADREAD datasets, which include RNA-seq expression matrices, clinical data, and masked annotated somatic mutations, were acquired from the UCSC repository. Additionally, the GSE39582 dataset, comprising 566 colorectal cancer (CRC) samples, and the GSE17538 dataset, containing 238 CRC samples with corresponding clinical information, were downloaded from the GEO database. Neddylation-related genes (NRGs) and gene sets from five pathways (Biocarta, Hallmark, KEGG, Reactome, and Wiki Pathways) were compiled from the Molecular Signatures Database (MsigDB, https://www.gsea-msigdb.org/gsea/msigdb/).
Single-Cell RNA-Seq Analysis
Two single-cell transcriptomic datasets, GSE166555 and GSE139555, were integrated for subsequent analysis utilizing the “Seurat” package,29 which is specifically designed for single-cell RNA sequencing analysis. To address batch effects, regularized negative binomial regression was employed. Non-linear dimensionality reduction was performed using Uniform Manifold Approximation and Projection (UMAP). Cluster-specific biomarkers were identified using the “COSG” package, while differential gene expression for each cell type was determined via the “FindAllMarkers” function. Furthermore, hallmark gene enrichment scores were evaluated using single-sample Gene Set Enrichment Analysis (ssGSEA) from the “GSVA” package, employing hallmark gene sets sourced from the Molecular Signatures Database (MsigDB).
Neddylation Patterns Identified via Unsupervised Clustering
We performed a univariate Cox regression analysis to identify NRGs that significantly impacted overall survival (OS) in CRC (p<0.01), utilizing the “survival”, “survminer”, and “ggplot2” packages. Subsequently, we employed the “ConsensusClusterPlus” package30 to perform unsupervised consensus clustering on the expression profiles of prognostic NRGs. The optimal number of clusters was determined based on the cumulative distribution curve and K-means clustering. The robustness of the clustering results was further validated through principal component analysis (PCA).
Signature of Neddylation-Related Genes (NRGS) Obtained Through Integrative Machine Learning methods
Based on the Univariate Cox regression analysis with a p-value threshold of <0.05, we identified 52 NRGs with significant prognostic value. To develop a robust prognostic NRG signature (NRGS), these biomarkers were subjected to an integrative machine learning analysis. This procedure included various algorithms such as random survival forests (RSF), least absolute shrinkage and selection operator (LASSO), gradient boosting machine (GBM), survival support vector machine (Survival-SVM), supervised principal components (SuperPC), ridge regression, partial least squares regression for Cox (plsRcox), CoxBoost, Stepwise Cox, and elastic network (Enet). Notably, RSF, LASSO, CoxBoost, and Stepwise Cox have capabilities for dimensionality reduction and variable selection, and they were combined with other algorithms to form multiple machine learning algorithm combinations. Similar methodologies have been employed in previous studies.31 The TCGA-COADREAD dataset was utilized as the training dataset, while the GSE17538 and GSE39582 datasets served as external validation datasets. The concordance index (C-index) for each model across all validation datasets was computed, and the model exhibiting the highest average C-index was identified as the optimal model. Following the identification of the optimal model, the median risk score derived from the training dataset was employed as the threshold to stratify patients in both the training and validation datasets into high-risk and low-risk groups. Kaplan-Meier survival analysis and the Log rank test were conducted on these two groups using the “survival” and “survminer” R packages. Receiver operating characteristic (ROC) curves were generated to assess the predictive accuracy of the model.
Gene Set Variation Analysis (GSVA)
The hallmark gene sets associated with neddylation and related pathways were obtained from the MsigDB. The enrichment scores of these hallmark genes were assessed using single-sample Gene Set Enrichment Analysis (ssGSEA) via the R package “GSVA.”32
Integrated Omics Analysis
The “maftools” package33 was utilized to characterize somatic mutations in genes of CRC patients.
Functional Enrichment Analysis
To explore the potential biological functions of differential expressed genes, Gene Ontology (GO) functions, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathway enrichment analyses were performed using GSEA with R package “clusterProfiler”.34
Immune Infiltration Analysis
The single-sample GSEA (ssGSEA) method was employed with the R package “GSVA” to evaluate immune cell scores. Furthermore, the R package “IOBR”35 which integrates eight algorithms (MCPcounter, EPIC, xCell, CIBERSORT, IPS, quanTIseq, ESTIMATE, and TIMER), was used to estimate the abundance of immune cells across different risk groups. Moreover, we analyzed hub gene expression in immune and non-immune cells within CRC single-cell datasets using the TISCH database.
Immunotherapeutic Response Prediction and Drug Sensitivity Analysis
To assess the role of NRGS in predicting the benefits of immunotherapy, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) score, accessible via the TIDE website (http://tide.dfci.harvard.edu). Subsequently, we investigated drug susceptibility across the two NRGS-defined risk groups by calculating the half-maximal inhibitory concentration (IC50) values using the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) with the “oncoPredict” R package.36
Statistical Analysis
Data processing, statistical analyses, and visualization were performed using R software (version 4.2.2). The Pearson correlation coefficient was utilized to evaluate the correlation between continuous variables. Chi-square tests were applied for comparisons of categorical variables, while Wilcoxon rank-sum tests or T-tests were used for continuous variables. The “survminer” package was employed to determine optimal cutoff values. Cox regression and Kaplan-Meier analyses were conducted using the “survival” package. Figures were generated with the “ggplot2” package. A p-value of less than 0.05 was considered statistically significant (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).
Results
The Allocation of Various Cell Subsets in CRC According to the scRNA-Seq Atlas
We utilized the scRNA-seq datasets (GSE166555 and GSE139555) and selected 14 samples from CRC patients for further investigation. Rigorous quality control procedures were implemented, resulting in the acquisition of 43,725 cells from the predetermined samples. Through UMAP analysis and hierarchical clustering, we identified 16 distinct cell subgroups, including B cells, CD4+ T cells, CD8+ T cells, endothelial cells, epithelial cells, fibroblasts, dendritic cells (DC), malignant cells, mast cells, monocytes/macrophages (Mono/Macro), myofibroblasts, natural killer (NK) cells, plasma cells, Th17 cells, proliferating T cells (Tprolif), and regulatory T cells (Treg) (Figure 1A). Subsequently, we examined the expression patterns of marker genes across these 16 cell types and presented the top three genes in a heat map (Figure 1B). Using the FindAllMarkers function, we identified differentially expressed genes for each of the 16 cell types. Figure 1C illustrates the top five up-regulated and down-regulated genes for each cell type. The GSVA R package was employed to compute the gene set variation analysis scores for 50 hallmark pathways across the cell types, revealing significant enrichment in endothelial, epithelial, fibroblasts, myofibroblasts, DC, Mono/Macro, and malignant cells (Figure 1D).
|
Figure 1 Different cell clustering based on scRNA-seq data of CRC and further analysis. (A) Cluster annotation and 16 cell types identification by means of UMAP. (B) Heat map of marker genes for the 16 cell types. (C) Differential expression genes for the 16 cell types. (D) Heat map of the scores of HALLMARK pathways for the 16 cell types. scRNA-seq, single cell RNA sequencing; UMAP, Unified Flowform Approximation and Projection.
|
The NRGs Status Across Various Cell Subsets as Determined by the scRNA-Seq Atlas
Utilizing the MSigDB database, a total of 250 neddylation-related genes (NRGs) were identified. Figure 2A illustrates the expression patterns of these 250 NRGs across 16 distinct cell types. Gene Set Variation Analysis (GSVA) was employed to compute the neddylation gene set scores for each cell type, categorizing them into high-score and low-score groups. The results were visualized using Uniform Manifold Approximation and Projection (UMAP) plots (Figure 2B). Notably, cell types such as Mono.Macro, DC, NK, malignant, endothelial, fibroblasts, myofibroblasts, and Tprolif exhibited relatively higher neddylation scores (Figure 2C). Figure 2D presents the composition of cell types within the two groups and across individual samples. Significant differences in cell type composition were observed between the high-score and low-score groups. The high-score group demonstrated increased proportions of malignant, Mono.Macro, fibroblasts, endothelial, myofibroblasts, DC, NK, Tprolif, and Treg cells, whereas the low-score group contained relatively higher proportions of CD4 T, B, and CD8 T cells. These findings were further elucidated through linear plots and UMAP plots, offering a comprehensive visualization of cell type distribution (Figure 2E and F).
 |
Figure 2 The status of NRG in different cell subsets based on scRNA-seq atlas. (A) The expression patterns of 250 NRG in 16 cell types. UMAP plots (B) and Box plots (C) showed differential NRG scores in 16 cell types. (D) Box plots showed the composition of the 16 cell types in high and low NRG score groups and in each sample. Linear plots (E) and UMAP plots (F) providing a comprehensive visualization of the 16 cell types distribution in two NRG score groups.
|
The HALLMARK Pathways Between High-and Low-Score Groups Based on scRNA-Seq Atlas
We subsequently conducted an analysis of the differential enrichment of hallmark pathways between groups with high and low scores. The findings indicated that nearly all pathways were significantly enriched in the high-score group, including the MYC targets pathway, MTORC1 signaling pathway, TGF-beta signaling pathway, KRAS signaling pathway, PI3K-AKT-MTOR signaling pathway, DNA repair, and oxidative phosphorylation pathways (Figure 3A). Furthermore, we examined the correlation between the neddylation score and hallmark pathways across various cell types. The results demonstrated that the neddylation score was positively correlated with almost all pathways in 16 cell types, particularly in endothelial, epithelial, fibroblast, dendritic, malignant, mast, monocyte/macrophage, myofibroblast, plasma, and proliferating T cells (Figure 3B). The most significant pathways included the MYC targets pathway, MTORC1 signaling pathway, PI3K-AKT-MTOR signaling pathway, DNA repair, and oxidative phosphorylation pathways.
 |
Figure 3 The HALLMARK pathways enrichment analysis in two NRG risk score groups based on scRNA-seq atlas. (A) The HALLMARK pathway enrichment between high and low NRG score groups. (B) Correlation of NRG score with HALLMARK pathways in 16 cell types. * p<0.05, ** p<0.01, *** p<0.001, **** p<0.0001.
|
Discovery of Various Neddylation Patterns in CRC Using Bulk-RNA Sequencing Data
Subsequently, we conducted a univariate Cox regression analysis on 250 NRGs in CRC samples, which led to the identification of 27 NRGs with significant prognostic value (p<0.01) that exhibited inter-correlation (Figure 4A). Utilizing a consensus clustering algorithm, optimal clustering was achieved at k=2 (Figure 4B), allowing us to categorize the CRC samples into two distinct neddylation patterns, designated as Cluster A (n=701) and Cluster B (n=455). Notably, patients in Cluster B demonstrated poorer overall survival compared to those in Cluster A (Figure 4C). Furthermore, there were significant differences in the expression levels of all prognostic NRGs, with the exception of COPS7A, between the two clusters (Figure 4D). Specifically, Cluster A exhibited higher expression levels of CISH, CCNF, BIRC5, PSME3, PSMB8, FBXW9, PSMA5, BRCA1, PSMA7, PSME1, DCAF7, PSMD3, PSMB10, and PSMB2, whereas Cluster B showed elevated expression of SOCS5, KLHL20, CCDC8, FBXL7, WSB1, FEM1B, COPS8, FBXO30, RNF7, COPS2, FBXO15, and FBXO32. Subsequently, we generated a heat map to visualize the clinicopathological factors and the expression of the 27 prognostic NRGs across the two clusters in CRC patients (Figure 4E). These findings support the hypothesis that distinct neddylation patterns are indicative of varying prognoses and clinical characteristics.
 |
Figure 4 Identification of two neddylation patterns of CRC based on bulk-RNA sequencing data. (A) The correlation network of 27 NRG with significant prognostic value in CRC identified by univariate Cox analysis with p-value <0.01 as the cut off. In the right half, purple denotes a prognostic risk factor, and green signifies a protective factor. Circle size reflects the P value magnitude. A line indicates a correlation between related genes with p < 0.05. (B) The optimal number of clusters based on k=2. (C) Kaplan-Meier analysis of prognosis between the two patterns. (D and E) Differences in clinical characteristics and prognostic NRG expression between the two patterns. * p<0.05, *** p<0.001, ns, not significant.
|
The Enrichment of Pathways Between the Two Neddylation Patterns
We subsequently acquired five pathway gene sets from the MSigDB database, specifically Biocarta, Hallmark, KEGG, Reactome, and WikiPathways. The R package “GSVA” was employed to score these pathways, and the R package “pheatmap” was utilized to generate a heatmap for comparing the two groups. Within the Biocarta pathway, Cluster A exhibited enrichment in pathways related to MCM, P27, G2, cell cycle, mitochondrial function, P53, MHC, Fas, and TNFR1. Conversely, Cluster B demonstrated enrichment in pathways associated with NK cells, CD40, TCRA, IL-6, ERK5, BCR, IL-10, IL-4, IL-3, MET, TCR, CXCR4, IL1R, IGF1, IL-7, PTEN, ECM, ALK, B lymphocytes, CCR5, TGF-β, LYM, and monocytes. Regarding the Hallmark pathway, Cluster A was enriched in MYC targets, G2M checkpoint, oxidative phosphorylation, DNA repair, and MTORC1 signaling. In contrast, Cluster B was enriched in IL6-JAK-STAT3 signaling, apoptosis, IL2-STAT5 signaling, NOTCH signaling, hypoxia, inflammatory response, TNFA signaling via NFKB, KRAS signaling, TGF-β signaling, angiogenesis, and EMT. In the context of KEGG pathways, Cluster A demonstrated significant enrichment in pathways related to DNA replication, mismatch repair, base excision repair, the cell cycle, and oxidative phosphorylation. Conversely, Cluster B exhibited enrichment in pathways associated with B cell receptor signaling, chemokine signaling, cytokine-cytokine receptor interactions, MAPK signaling, adherens junctions, TGF-β signaling, focal adhesion, and ECM receptor interactions. Regarding Reactome pathways, Cluster A was enriched in processes such as DNA replication, mismatch repair, G2-M checkpoints, base excision repair, cell cycle checkpoints, G2-M DNA damage checkpoints, stabilization of P53, DNA repair, cellular response to hypoxia, regulation of TP53 activity, apoptosis, and MTORC1-mediated signaling. In contrast, Cluster B was enriched in pathways involving IL-3, IL-5, and GM-CSF signaling, IL-27 signaling, IL-21 signaling, PI3K-AKT signaling in cancer, RET signaling, IL-15 signaling, AKT signaling in cancer, and the IL-6 signaling pathway. In the WikiPathway analysis, Cluster A exhibited significant enrichment in pathways related to DNA replication, DNA mismatch repair, base excision repair, metabolic reprogramming in colon cancer, the cell cycle, and ATM signaling. Conversely, Cluster B demonstrated enrichment in pathways associated with TGF-β signaling, B cell receptor signaling, T cell receptor signaling, NOTCH signaling, chemokine signaling, IL-3 signaling, Wnt signaling, IL-4 signaling, MAPK signaling, RAS signaling, IL-10 anti-inflammatory signaling, PI3K-AKT signaling, focal adhesion, TGF-β receptor signaling, CCL18 signaling, and angiogenesis. These findings indicate distinct immune activities and proliferative statuses between the two risk groups, potentially explaining the observed differences in survival rates (Figure 5).
 |
Figure 5 Pathways enrichment analysis in two neddylation patterns. The pathways enrichment of Biocarta (A), Hallmark (B), KEGG (C), Reactome (D) and wiki pathway (E) between the two patterns.
|
Analyses of Somatic Mutation in Two Neddylation Patterns
We subsequently examined the differences in genomic mutations associated with the two neddylation patterns. The somatic mutation landscape of each CRC sample within the cluster A and B groups was visualized using a waterfall diagram. It was observed that the top ten mutant genes in the two subgroups were largely consistent (Figure 6A and B). In the cluster A group, the genes with the highest mutation frequencies were APC (61%), TP53 (57%), TTN (50%), KRAS (46%), and MUC16 (28%). Similarly, in the cluster B group, the most frequently mutated genes were APC (60%), TP53 (55%), KRAS (42%), TTN (40%), and SYNE1 (27%). Regarding the classification of variations, missense mutations, nonsense mutations, and multi-hit mutations were the three most prevalent types across all mutations. We further analyzed the differential mutant genes between the two groups. Figure 6C illustrates the top ten genes with higher mutation frequencies in the cluster A group compared to the cluster B group, namely RIMS1, SPEN, FAT1, FSIP2, YTHDC1, ZNF516, MXRA5, ASTN2, PCNT, and PCDHGA6. Additionally, tumor mutational burden (TMB) was compared between the two groups, and no significant difference was observed (data not shown).
 |
Figure 6 Somatic mutation analysis in two neddylation patterns. (A) The top 10 mutated genes in the cluster A group. (B) The top 10 mutated genes in the cluster B group. (C) Differences in mutated genes between the two patterns. ** p<0.01, *** p<0.001.
|
The Immune Landscapes Differ Between the Two Neddylation Patterns
Principal Component Analysis (PCA) demonstrated that the various subtypes were heterogeneous and potentially indicative of differing levels of neddylation modifications (Figure 7A). Furthermore, these findings were corroborated using single-sample Gene Set Enrichment Analysis (ssGSEA), revealing that the cluster A group typically exhibited a lower abundance of immune cells, including activated B cells, activated dendritic cells, macrophages, mast cells, Myeloid-derived suppressor cells (MDSCs), NK cells, NKT cells, regulatory T cells, and type 1 T helper cells (Figure 7B). Subsequently, eight algorithms were employed to investigate the variations in the tumor microenvironment (TME) status across the two neddylation patterns. The analysis indicated that the cluster A group had lower immune, stromal, and ESTIMATE scores, but a higher tumor purity score. This group also showed a relatively low abundance of infiltrated immune cell types (Figure 7C). In summary, both analytical approaches suggested that cluster A exhibited reduced levels of immune infiltration, which may contribute to tumor immune evasion.
 |
Figure 7 Different immune landscapes of two neddylation patterns. (A) PCA in two patterns. (B) Abundance of 23 immune cells based on ssGSEA algorithm in two patterns. (C) Eight immune infiltration algorithms exhibit different numbers of immune cells between the two patterns. * p<0.05, ** p<0.01, *** p<0.001, ns, not significant.
|
Development and Verification of Prognostic Signatures Associated with Neddylation-Related Genes (NRGS)
Utilizing a Univariate Cox regression analysis with a p-value threshold of <0.05, we identified 52 NRGs with significant prognostic implications, comprising 24 risk-associated genes and 28 protective genes. These 52 potential prognostic biomarkers were subsequently subjected to an integrative machine-learning procedure involving 10 previously mentioned methods. We computed the average concordance index (C-index) for each algorithm across all cohorts and selected the StepCox[both]+RSF algorithm, which exhibited the highest average C-index of 0.71, as the optimal model (Figure 8A). Subsequently, we calculated the risk score for each sample within the cohort based on the expression profiles of the 26 NRGs included in the StepCox[both]+RSF model, namely PLAB2, DCAF7, PSMA5, ASB9, COPS7A, FBXO15, COPS8, FBXL16, RNF7, CCDC8, PSMD12, FBXL19, BIRC5, BTRC, PSME1, TRIM40, FBXL7, DCAF5, CCNF, CUL4B, NAE1, KLHL21, SOCS3, DCUN1D3, FBXO32, and OBSL1 (Figure 8B). Kaplan-Meier survival analysis for overall survival (OS) demonstrated that the high-risk score group exhibited significantly poorer survival outcomes in the training cohort (TCGA-CRC) (Figure 8C). Importantly, the area under the receiver operating characteristic (ROC) curve (AUC) for 1-year, 2-year, and 3-year OS was 0.981, 0.987, and 0.987, respectively, in the training cohort (Figure 8D). Moreover, we showed the changes in AUC at different time points in the training cohort and found that the AUC is all close to 1 (Figure 8E). High-score group also had poorer overall survival in the testing cohort (GSE17538 and GSE39582) and all combined cohorts (Figure 8F–H).
 |
Figure 8 A prognostic NRGS developed by machine learning analysis. (A) The C-index of 107 kinds of prognostic models constructed by 10 machine learning algorithms in training and testing cohort. (B) Error rate curve of random forest tree model. (C) Kaplan-Meier analysis of different NRGS score in training cohort (TCGA-COADREAD). (D) The area under the ROC curve (AUC) for the 1-year, 2-year, and 3-year OS in the training cohort. (E) The changes in AUC value of different time points. Kaplan-Meier analysis of different NRGS score in testing cohort (GSE17538) (F), (GSE39582) (G) and in merge cohort (H).
|
We also investigated the differences in genomic mutations between the high- and low- risk groups. TP53 (60%), APC (53%), KRAS (52%), TTN (45%) and SYNE1 (29%) were the top 5 genes with the highest mutation frequencies in the high-risk group (Figure S1A), while APC (61%), TP53 (56%), TTN (46%), KRAS (43%) and SYNE1 (27%) in the low risk groups (Figure S1B). Missense mutation, nonsense mutation and multi hit are also the top 3 across all mutation types. Figure S1C showed the top 10 higher mutant genes (ENPEP, DLCK1, DCAF12, GAREML, USP12, SUN2, CD74, LYNX1, ADORA3 and GAGNG3) in high-risk group compared with low-risk group.
We utilized the TISCH database to investigate the distribution of the top ten genes (ASB9, CCDC8, COPS7A, COPS8, DCAF7, FBXL16, FBXO15, PALB2, PSMA5, and RNF7) associated with the RSF model at the single-cell level (Figure S2). Our analysis revealed that the distributions of COPS7A, COPS8, DCAF7, PSMA5, and RNF7 were widespread across most detected cell types, particularly in Tprolif cells.
Correlation Between Clinical Features and NRGS
Furthermore, significant differences were observed in the distribution of high and low-score groups concerning recurrence or metastasis, stage, and tumor type (Figure 9A–D). Subsequently, we categorized all patients based on various clinical characteristics and found that the risk score significantly increased with recurrence or metastasis, advanced stage, and COAD (Figure 9A–D). Notably, as illustrated in Figure 9E and F, both univariate and multivariate Cox regression analyses demonstrated that the NRGS-based risk score served as an independent risk factor for the overall survival rate of CRC.
 |
Figure 9 Correlation between Clinical Features and NRGS. The differences in risk scores across clinical subgroups, including gender (A), recurrence or metastasis (B), stage (C) and tumor type (D). Univariate (E) and multivariate Cox regression analysis (F) illustrated that the NRGS could be used as an independent prognostic factor for CRC patients (P< 0.001).
|
Functional Enrichment Analysis of NRGS
GSEA was employed to further elucidate the biological functions and signaling pathways associated with NRGS. Initially, we examined the genes correlated with the NRGS-based risk score in CRC. Figure 10A and B illustrate the top 50 genes demonstrating positive and negative correlations with the risk score, respectively. Regarding GO terms, genes related to the NRGS risk score were enriched in processes such as extracellular matrix organization, extracellular structure organization, collagen metabolic process, blood vessel development, angiogenesis, negative regulation of tumor necrosis factor production, cell-substrate adhesion, cell-matrix adhesion, and negative regulation of tumor necrosis factor superfamily cytokine production (Figure 10C). In terms of the KEGG pathways, genes associated with the NRGS risk score were enriched in ECM-receptor interaction, phagosome, focal adhesion, cell adhesion molecules, leukocyte transendothelial migration, protein digestion and absorption, calcium signaling pathway, and proteoglycans in cancer (Figure 10D). For Reactome pathways, NRGS risk score-related genes showed enrichment in extracellular matrix organization, amplification of signal from the kinetochores, mitotic spindle checkpoint, G2/M checkpoints, and DNA strand elongation (Figure 10E).
 |
Figure 10 Functional Enrichment Analysis of NRGS. (A) The top 50 genes positively correlated with NRGS score. (B) The top 50 genes negatively correlated with NRGS score. (C) GO functional enrichment analysis of NRGS score related genes. (D) KEGG pathway enrichment analysis of NRGS score related genes. (E) Reactome pathway enrichment analysis of NRGS score related genes.
|
The TME Status Differ Between the Two NRGS Risk Groups
We employed eight algorithms to investigate the variations in the TME status between high-risk and low-risk groups. The findings from the ESTIMATE algorithm indicated that stromal, immune, and ESTIMATE scores were significantly elevated in the high-risk group compared to the low-risk group, whereas tumor purity was notably higher in the low-risk group. Additionally, cell populations such as activated dendritic cells (aDC), cancer-associated fibroblasts (CAFs), central memory CD4+ T cells (CD4-Tcm), conventional dendritic cells (cDC), cytotoxic lymphocytes, dendritic cells (DC), endothelial cells, immature dendritic cells (IDC), macrophages, M0 macrophages, M2 macrophages, major histocompatibility complex molecules (MHC), monocytes, mesenchymal stem cells (MSC), myeloid dendritic cells, neutrophils, and CD8+ T cells were more prevalent in the high-risk group. Conversely, cell types such as CD4+ memory T cells, CD4+ naive T cells, immunosuppressive cells (SC), immune checkpoint pathways (IPS), resting natural killer (NK) cells, natural killer T (NKT) cells, plasmacytoid dendritic cells (pDC), T helper 1 (Th1) cells, T helper 2 (Th2) cells, and azurophilic granules (AZ) were more predominant in the low-risk group (Figure 11A).
 |
Figure 11 Correlation between immune microenvironment and NRGS in CRC. (A) The correlation between NRGS and the immune cell infiltration is based on eight algorithms. (B) The expression of chemokines and receptors, interieuknes and receptors, interferons and receptors, and other cytokines in CRC patients with different NRGS scores. * p<0.05, ** p<0.01, *** p<0.001.
|
Subsequently, we conducted an analysis of the differential expression levels of chemokines and their receptors, interleukins and their receptors, interferons and their receptors, as well as other cytokines between high- and low-risk groups. Regarding chemokines and their receptors, the expression levels of CCL5, CCL8, CCL17, CCL18, CCL21, CCL22, CCR1, CCR2, CCR5, CCR10, CXCL12, CXCL16, and CXCR4 were elevated in the high-risk group, whereas CCL20, CCL28, CCR6, and CXCL11 exhibited higher expression in the low-risk group. In terms of interleukins and their receptors, IL1R1, IL1R2, IL4R, IL6, IL10, IL10RA, IL10RB, IL11, IL17B, IL17D, and IL21R showed increased expression in the high-risk group, while IL1A, IL5, IL12A, and IL26 were more highly expressed in the low-risk group. For interferons and their receptors, IFNAR2 and IFNGR2 were more highly expressed in the high-risk group. Among other cytokines, CSF1, CSF2RB, PDGFC, PDGFD, PDGFRA, PDGFRB, TGFB3, TGFBR1, TGFBR2, TNF, VEGFB, VEGFC, and EPOR exhibited higher expression levels in the high-risk group, with EGF being the only cytokine with elevated expression in the low-risk group (Figure 11B).
Analysis of Predicted Responses to Immunotherapy of NRGS
We assessed the Tumour Immune Dysfunction and Exclusion (TIDE) scores of CRC patients using the TIDE database. Our analysis revealed that the high-risk group exhibited significantly higher TIDE scores (Figure 12A). Bar plots depicting the distribution of responders and non-responders within both high- and low-risk groups indicated that 75% of individuals in the high-risk group were non-responders, while 25% were responders. In contrast, the low-risk group comprised 51% non-responders and 49% responders (Figure 12B). Furthermore, the high-risk group demonstrated an increased propensity for tumour immune exclusion (Figure 12C) and dysfunction (Figure 12D).
 |
Figure 12 Prediction of immunotherapy effects of NRGS in CRC. (A) The TIDE score in high- and low-risk groups. (B) The immunotherapy response rate in high- and low-risk groups. The T-cell exclusion score (C) and T-cell dysfunction score (D) in high- and low-risk groups.
|
Drug Sensitivity of NRGS
Subsequently, the oncoPredict algorithm was employed to estimate the IC50 values, facilitating the prediction of differential chemotherapy responses between high-risk and low-risk patient groups. A significant variance in the sensitivity to numerous chemotherapeutic agents was observed between these groups. Specifically, the results indicated that patients in the high-risk group exhibited lower sensitivity to oxaliplatin (p<0.0001), 5-Fluorouracil (p<0.01), KRAS G12C inhibitor (p<0.0001), sorafenib (p<0.0001), trametinib (p<0.05), dabrafenib (p<0.01), afatinib (p<0.0001), cisplatin (p<0.05), and crizotinib (p<0.001) compared to those in the low-risk group, as illustrated in Figure 13.
 |
Figure 13 Prediction of drug sensitivity of NRGS in CRC. The Y-axis showed the IC50 value, which was negatively correlated with drug sensitivity. *p<0.05, ** p<0.01, ****p < 0.0001.
|
Discussion
Patients with advanced metastatic colorectal cancer (mCRC) frequently miss the opportunity for surgical intervention, rendering pharmacological treatment the primary therapeutic strategy. In recent years, immune checkpoint inhibitors (ICIs) have demonstrated remarkable and sustained survival benefits in a limited subset of patients characterized by high immunogenicity, specifically those with microsatellite instability-high (MSI-H) or mismatch repair-deficient (dMMR) tumors.3 Conversely, the majority of colorectal cancer (CRC) patients exhibit low immunogenicity, such as microsatellite instability-low (MSI-L) or proficient mismatch repair (pMMR), and consequently show resistance to ICIs. Furthermore, numerous studies have sought to identify several distinct CRC subtypes to predict the prognosis and efficacy of ICIs treatment, including those classified by the Consensual Molecular Subtype (CMS) classification, as well as mutations in DNA polymerase D1 (POLD1) or DNA polymerase E (POLE), and Lynch syndrome (LS), among others.4 Therefore, there remains a critical and unmet need to identify molecular subtypes or biomarkers that can enhance the selection of patients likely to respond to specific treatments. Neddylation, a post-translational modification, is frequently hyperactive in cancer and significantly impacts the tumor microenvironment (TME) by modulating immune cell activities, positioning it as a promising target for cancer therapy.19,20 However, systematic investigations into neddylation-related genes (NRGs) in CRC are currently lacking. Therefore, we conducted a systematic analysis of NRGs in CRC to elucidate their potential roles in tumor progression, prognosis, the TME and treatment responses. This was achieved through the characterization of the single-cell landscape and the integration of bulk transcriptome data.
Initially, we integrated single-cell data from the GSE166555 and GSE139555 cohorts to analyze cellular heterogeneity in CRC. Our analysis indicated that the group with a high neddylation score exhibited increased proportions of malignant cells, monocytes/macrophages, fibroblasts, endothelial cells, myofibroblasts, DC cells, NK cells, and Tprolif cells, all of which are essential components of the tumor microenvironment (TME). These TME components are instrumental in facilitating in tumor progression, and targeting these cells could be crucial in influencing tumor outcomes. Notably, Tprolif cells exhibited the highest neddylation scores and were characterized by elevated expression levels of immunological exhaustion markers, including PDCD1, HAVCR2, CTLA4, LAG3, and TIGIT.37 Importantly, an increasing body of research underscores the significant role of neddylation in regulating the TME.19,20 These findings suggest that neddylation may significantly influence tumor-stroma interactions in CRC.38
In a subsequent analysis, we identified 27 neddylation-related genes (NRGs) with significant prognostic value (p<0.01) through univariate Cox regression analysis. Employing a consensus clustering algorithm, we further delineated two distinct neddylation patterns, termed Cluster A and Cluster B. Notably, differential expression was observed in all prognostic NRGs, with the exception of COPS7A, between these two clusters. Cluster B was associated with a poorer prognosis. Single-sample Gene Set Enrichment Analysis (ssGSEA) indicated that Cluster B exhibited a higher abundance of immune cells, including B cells, T cells, and natural killer (NK) cells. However, it also exhibited an increased presence of cells functioning as immunosuppressive regulators in cancer, such as myeloid-derived suppressor cells (MDSCs),39 regulatory T cells (Tregs),40 natural killer T (NKT) cells,41 and tumor-associated macrophages (TAMs).42–44 Interestingly, Treg cells have been observed to enhance TAM activity, thereby inducing a suppressive TME through the inhibition of CD8+ T cells.45 This cellular composition may contribute to a more aggressive cancer phenotype and a less favorable prognosis for patients. Additionally, Cluster B demonstrated a significant enrichment in several signaling pathways, including TGF-β, IL6-JAK-STAT3, PI3K-AKT, NOTCH, MAPK, RAS, and angiogenesis pathways. Extensive research has established that the aberrant activation of these pathways plays a pivotal role in the proliferation, migration, and invasion of CRC, thus altering the TME and contributing to a poorer prognosis.46–51 Notably, neddylation is crucial in modulating the TME by influencing various signaling pathways, such as TGF-β, PI3K-AKT-mTOR, and EGFR pathways.52 These findings suggest that neddylation may facilitate tumor growth and metastasis by modulating the TME through the regulation of numerous signaling pathways in CRC.
To improve the accuracy of risk stratification in CRC patients, a neddylation-related gene signature (NRGS) was developed. This signature demonstrated high predictive accuracy, with patients exhibiting higher risk scores experiencing significantly poorer outcomes across the training, validation, and combined cohorts. The area under the curve (AUC) for 1-year, 2-year, and 3-year overall survival (OS) in the training cohort was 0.979, 0.989, and 0.996, respectively. The NRGS-based risk score was associated with an increased likelihood of recurrence or metastasis, more advanced disease stages, and functioned as an independent prognostic factor for overall survival (OS) in CRC. Genetic mutation analysis indicated that patients with elevated risk scores exhibited a higher mutation frequency compared to those with lower risk scores, suggesting that an increase in genetic mutations may lead to cellular physiological dysfunction and promote tumor metastasis. Notably, the group characterized by elevated NRGS demonstrated a higher prevalence of immunosuppressive cells, including myeloid-derived suppressor cells (MDSCs), mesenchymal stem cells (MSCs), and M2-tumor-associated macrophages (TAMs). Furthermore, this group exhibited an upregulation of chemokines and their corresponding receptors. The intricate interactions among dysregulated cytokines, chemokines, growth factors, and matrix-remodeling enzymes play a critical role in to the pathogenesis of CRC and elicit systemic responses that influence disease outcomes.53 The upregulation of CCR/CCL5 expression has been correlated with the infiltration of immune cells that are associated with poor prognosis, including regulatory T cells (Tregs), M1 and M2 macrophages, myeloid-derived suppressor cells, and cancer-associated fibroblasts.54 Additionally, increased CCR/CCL5 expression levels are associated with a wide array of immunosuppressive proteins, such as PD-1, PD-L1, PD-L2, CTLA4, CD80, CD86, TIM3, IDO1, LAG3, and IFN-γ, suggesting potential mechanisms by which CRC circumvents anti-cancer immune responses.54 Macrophage-derived CCL5 facilitates the immune evasion of colorectal cancer cells via the p65/STAT3-CSN5-PD-L1 pathway.55 Moreover, the neddylation pathway enhances CCL2 transactivation and tumor-associated macrophage (TAM) infiltration,56 which are critical for tumor immunosuppression, thereby creating a microenvironment that supports to cancer development and progression.
Significantly, the cohort exhibiting elevated NRGS demonstrated a low IPS score, a low TIDE score, and increased scores for tumor immune exclusion and dysfunction. The IPS has been identified as a superior predictor of response to anti-CTLA4 and anti-PD1 antibodies, with a low IPS score suggesting a diminished response to immunotherapy.57 Conversely, a low TIDE score indicates a reduced likelihood of immune escape and an enhanced response to immunotherapy.58,59 Furthermore, the high NRGS group exhibited resistance to standard colorectal cancer therapies, including oxaliplatin, 5-Fluorouracil, KRAS G12C inhibitors, and sorafenib. Therefore, this signature could help clinicians identify patient subgroups for tailored immunotherapy and chemotherapy, but larger patient groups are needed for validation.
Limitations of the Study
This study represents the inaugural effort to identify prognostically significant neddylation-related genes in colorectal cancer (CRC) through the utilization of scRNA-seq and bulk RNA sequencing data. Nonetheless, several limitations are evident: firstly, the sample size for scRNA-seq was insufficient, and there was notable heterogeneity between scRNA-seq and bulk RNA-seq data. Secondly, although the Neddylation-Related Gene Signature (NRGS) demonstrated efficacy in predicting prognosis in CRC patients, further validation with larger and more diverse clinical samples is necessary to enhance the reliability of the model’s predictive capacity. Additionally, the study’s reliance on public databases without the inclusion of clinical samples, the absence of validation of gene expression in CRC tissues, and the lack of cellular or animal experiments to evaluate the genes’ impact on CRC behavior, constitute significant limitations. Consequently, further research is warranted to elucidate the roles and mechanisms of these genes in CRC progression. Despite these constraints, the identified markers provide valuable guidance for future investigations and hold potential for informing novel therapeutic targets and clinical strategies for CRC.
Conclusion
This study conducted a comprehensive investigation into the roles of neddylation-related genes (NRGs) in the progression of colorectal cancer (CRC), as well as their impact on responses to immunotherapy and chemotherapy. Our results indicated that patients with a high neddylation score exhibited a greater presence of malignant cells and diverse immune cell populations, alongside activation of critical pathways associated with tumor proliferation and immune evasion. We developed a neddylation-related gene signature (NRGS), which effectively characterizes the immunological landscape of CRC patients and provides a reliable prognostic tool. This signature can aid clinicians in identifying specific patient subgroups that may derive benefit from tailored immunotherapy and chemotherapy regimens. Future research should prioritize the validation of these findings in larger patient cohorts and investigate the therapeutic potential of targeting NRGs in the treatment of CRC.
Data Sharing Statement
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Ethics Approval and Consent to Participate
All relevant ethical regulations were followed the original study of the datasets and the authors of the source studies had also obtained informed consent from participants. Ethical approval for this study was exempted by the Fourth Affiliated Hospital of Guangzhou Medical University, as the data were obtained from public sources.
Author Contributions
All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.
Funding
This study were supported by Guangzhou Zengcheng District Science and Technology Innovation fund project (ZCKJ2019-008 and 2021ZCMS18), and Science and Technology Program of Guangzhou, China (2024A03J0948).
Disclosure
The authors declare that they have no competing interests in this work.
References
1. Siegel RL, Wagle NS, Cercek A, Smith RA, Jemal A. Colorectal cancer statistics, 2023. CA: a Cancer Journal for Clinicians. 2023;73(3):233–254. doi:10.3322/caac.21772
2. Ganesh K, Stadler ZK, Cercek A, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol. 2019;16(6):361–375. doi:10.1038/s41575-019-0126-x
3. Ganesh K. Optimizing immunotherapy for colorectal cancer. Nature Reviews. Gastroenterology & Hepatology. 2022;19(2):93–94. doi:10.1038/s41575-021-00569-4
4. Zhou C, Cheng X, Tu S. Current status and future perspective of immune checkpoint inhibitors in colorectal cancer. CANCER LETT. 2021;521:119–129. doi:10.1016/j.canlet.2021.07.023
5. Zhang S, Yu Q, Li Z, Zhao Y, Sun Y. Protein neddylation and its role in health and diseases. Signal Transduct Target Ther. 2024;9(1):85. doi:10.1038/s41392-024-01800-9
6. Enchev RI, Schulman BA, Peter M. Protein neddylation: beyond cullin-RING ligases. Nat Rev Mol Cell Biol. 2015;16(1):30–44. doi:10.1038/nrm3919
7. Li L, Wang M, Yu G, et al. Overactivated neddylation pathway as a therapeutic target in lung cancer. J Natl Cancer Inst. 2014;106(6):u83. doi:10.1093/jnci/dju083
8. Jones TM, Carew JS, Bauman JE, Nawrocki ST. Targeting NEDDylation as a Novel Approach to Improve the Treatment of Head and Neck Cancer. Cancers (Basel). 2021;13(13):3250. doi:10.3390/cancers13133250
9. Norton JP, Augert A, Eastwood E, Basom R, Rudin CM, MacPherson D. Protein neddylation as a therapeutic target in pulmonary and extrapulmonary small cell carcinomas. GENE DEV. 2021;35(11–12):870–887. doi:10.1101/gad.348316.121
10. Kuo KL, Ho I-L, Shi C-S. MLN4924, a novel protein neddylation inhibitor, suppresses proliferation and migration of human urothelial carcinoma: in vitro and in vivo studies. CANCER LETT. 2015;363(2):127–136. doi:10.1016/j.canlet.2015.01.015
11. Jin Y, Zhang P, Wang Y. Neddylation Blockade Diminishes Hepatic Metastasis by Dampening Cancer Stem-Like Cells and Angiogenesis in Uveal Melanoma. CLIN CANCER RES. 2018;24(15):3741–3754. doi:10.1158/1078-0432.CCR-17-1703
12. Yao WT, Wu J-F, Yu G-Y. Suppression of tumor angiogenesis by targeting the protein neddylation pathway. CELL DEATH DIS. 2014;5(2):e1059. doi:10.1038/cddis.2014.21
13. Luo Z, Yu G, Lee HW. The Nedd8-activating enzyme inhibitor MLN4924 induces autophagy and apoptosis to suppress liver cancer cell growth. CANCER RES. 2012;72(13):3360–3371. doi:10.1158/0008-5472.CAN-12-0388
14. Yang Z, Zhang J, Lin X. Inhibition of neddylation modification by MLN4924 sensitizes hepatocellular carcinoma cells to sorafenib. ONCOL REP. 2019;41(6):3257–3269. doi:10.3892/or.2019.7098
15. Li H, Zhou W, Li L. Inhibition of Neddylation Modification Sensitizes Pancreatic Cancer Cells to Gemcitabine. NEOPLASIA. 2017;19(6):509–518. doi:10.1016/j.neo.2017.04.003
16. Shapiro DD, Zacharias NM, Tripathi DN. Neddylation inhibition sensitises renal medullary carcinoma tumours to platinum chemotherapy. Clin Transl Med. 2023;13(5):e1267. doi:10.1002/ctm2.1267
17. Wei D, Li H, Yu J. Radiosensitization of human pancreatic cancer cells by MLN4924, an investigational NEDD8-activating enzyme inhibitor. CANCER RES. 2012;72(1):282–293. doi:10.1158/0008-5472.CAN-11-2866
18. Vanderdys V, Allak A, Guessous F. The Neddylation Inhibitor Pevonedistat (MLN4924) Suppresses and Radiosensitizes Head and Neck Squamous Carcinoma Cells and Tumors. MOL CANCER THER. 2018;17(2):368–380. doi:10.1158/1535-7163.MCT-17-0083
19. Zhou L, Jiang Y, Luo Q, Li L, Jia L. Neddylation: a novel modulator of the tumor microenvironment. MOL CANCER. 2019;18(1):77. doi:10.1186/s12943-019-0979-1
20. Mao H, Lin X, Sun Y. Neddylation Regulation of Immune Responses. Research (Washington). 2023;6:283.
21. Zhou S, Zhao X, Yang Z. Neddylation inhibition upregulates PD‐L1 expression and enhances the efficacy of immune checkpoint blockade in glioblastoma. INT J CANCER. 2019;145(3):763–774. doi:10.1002/ijc.32379
22. Cui Y, Chen Z, Pan B. Neddylation pattern indicates tumor microenvironment characterization and predicts prognosis in lung adenocarcinoma. Front Cell Dev Biol. 2022;10:979262. doi:10.3389/fcell.2022.979262
23. Li Y, Niu J, Wang Y. Machine learning-based neddylation landscape indicates different prognosis and immune microenvironment in endometrial cancer. Frontiers in Oncology. 2023;13:1084523. doi:10.3389/fonc.2023.1084523
24. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. NAT METHODS. 2008;5(7):621–628. doi:10.1038/nmeth.1226
25. Tang F, Barbacioru C, Wang Y. mRNA-Seq whole-transcriptome analysis of a single cell. NAT METHODS. 2009;6(5):377–382. doi:10.1038/nmeth.1315
26. Tang Q, Tang L, Wang X. Comprehensive Analyses of Single-Cell and Bulk RNA Sequencing Data From M2 Macrophages to Elucidate the Immune Prognostic Signature in Patients with Gastric Cancer Peritoneal Metastasis. Immunotargets Ther. 2025;14:383–402. doi:10.2147/ITT.S506143
27. Tao Y, Li J, Pan J. Integration of scRNA-Seq and Bulk RNA-Seq Identifies Circadian Rhythm Disruption-Related Genes Associated with Prognosis and Drug Resistance in Colorectal Cancer Patients. Immunotargets Ther. 2025;14:475–489. doi:10.2147/ITT.S499806
28. Asada K, Takasawa K, Machino H. Single-Cell Analysis Using Machine Learning Techniques and Its Application to Medical Research. Biomedicines. 2021;9(11):1513. doi:10.3390/biomedicines9111513
29. Sun Y, Wu L, Zhong Y. Single-cell landscape of the ecosystem in early-relapse hepatocellular carcinoma. CELL. 2021;184(2):404–421. doi:10.1016/j.cell.2020.11.041
30. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. BIOINFORMATICS. 2010;26(12):1572–1573. doi:10.1093/bioinformatics/btq170
31. Liu Z, Liu L, Weng S. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. NAT COMMUN. 2022;13(1):816. doi:10.1038/s41467-022-28421-6
32. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC BIOINFORMATICS. 2013;14(7). doi:10.1186/1471-2105-14-7
33. Mayakonda A, Lin D, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. GENOME RES. 2018;28(11):1747–1756. doi:10.1101/gr.239244.118
34. Wu T, Hu E, Xu S. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3):100141. doi:10.1016/j.xinn.2021.100141
35. Zeng D, Ye Z, Shen R. IOBR: multi-Omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and Signatures. Front Immunol. 2021;12:687975. doi:10.3389/fimmu.2021.687975
36. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021;22(6). doi:10.1093/bib/bbab260
37. Yang C, Guo Y, Wu Z, Huang J, Xiang B. Comprehensive Analysis of Cuproptosis-Related Genes in Prognosis and Immune Infiltration of Hepatocellular Carcinoma Based on Bulk and Single-Cell RNA Sequencing Data. Cancers (Basel). 2022;14(22):5713. doi:10.3390/cancers14225713
38. Olaizola P, Lee-Law PY, Fernandez-Barrena MG. Targeting NAE1-mediated protein hyper-NEDDylation halts cholangiocarcinogenesis and impacts on tumor-stroma crosstalk in experimental models. J HEPATOL. 2022;77(1):177–190. doi:10.1016/j.jhep.2022.02.007
39. Li K, Shi H, Zhang B. Myeloid-derived suppressor cells as immunosuppressive regulators and therapeutic targets in cancer. Signal Transduction and Targeted Therapy. 2021;6(1):362. doi:10.1038/s41392-021-00670-9
40. Tanaka A, Sakaguchi S. Regulatory T cells in cancer immunotherapy. CELL RES. 2017;27(1):109–118. doi:10.1038/cr.2016.151
41. Seino K, Taniguchi M. Functionally distinct NKT cell subsets and subtypes. The Journal of Experimental Medicine. 2005;202(12):1623–1626. doi:10.1084/jem.20051600
42. Hong SM, Lee A-Y, Kim B-J. NAMPT-Driven M2 Polarization of Tumor-Associated Macrophages Leads to an Immunosuppressive Microenvironment in Colorectal Cancer. Adv Sci (Weinh). 2024;11(14):e2303177. doi:10.1002/advs.202303177
43. Huang C, Ou R, Chen X. Tumor cell-derived SPON2 promotes M2-polarized tumor-associated macrophage infiltration and cancer progression by activating PYK2 in CRC. J Exp Clin Cancer Res. 2021;40(1):304. doi:10.1186/s13046-021-02108-0
44. Han S, Wang W, Wang S. Tumor microenvironment remodeling and tumor therapy based on M2-like tumor associated macrophage-targeting nano-complexes. THERANOSTICS. 2021;11(6):2892–2916. doi:10.7150/thno.50928
45. Liu C, Chikina M, Deshpande R. Treg Cells Promote the SREBP1-Dependent Metabolic Fitness of Tumor-Promoting Macrophages via Repression of CD8(+) T Cell-Derived Interferon-gamma. IMMUNITY. 2019;51(2):381–397. doi:10.1016/j.immuni.2019.06.017
46. Prossomariti A, Piazzi G, Alquati C, Ricciardiello L. Are Wnt/beta-Catenin and PI3K/AKT/mTORC1 Distinct Pathways in Colorectal Cancer? Cell Mol Gastroenterol Hepatol. 2020;10(3):491–506. doi:10.1016/j.jcmgh.2020.04.007
47. Calon A, Espinet E, Palomo-Ponce S. Dependency of colorectal cancer on a TGF-beta-driven program in stromal cells for metastasis initiation. CANCER CELL. 2012;22(5):571–584. doi:10.1016/j.ccr.2012.08.013
48. Jackstadt R, van Hooff SR, Leach JD. Epithelial NOTCH Signaling Rewires the Tumor Microenvironment of Colorectal Cancer to Drive Poor-Prognosis Subtypes and Metastasis. CANCER CELL. 2019;36(3):319–336. doi:10.1016/j.ccell.2019.08.003
49. Pennel K, Hatthakarnkul P, Wood CS. JAK/STAT3 represents a therapeutic target for colorectal cancer patients with stromal-rich tumors. J Exp Clin Cancer Res. 2024;43(1):64. doi:10.1186/s13046-024-02958-4
50. Kotsiliti E. MAPK inhibition in BRAFV600E CRC. Nature Reviews. Gastroenterology & Hepatology. 2023;20(4):201. doi:10.1038/s41575-023-00756-5
51. Zhu G, Pei L, Xia H, Tang Q, Bi F. Role of oncogenic KRAS in the prognosis, diagnosis and treatment of colorectal cancer. MOL CANCER. 2021;20(1):143. doi:10.1186/s12943-021-01441-4
52. Liu D, Che X, Wu G. Deciphering the role of neddylation in tumor microenvironment modulation: common outcome of multiple signaling pathways. Biomark Res. 2024;12(5). doi:10.1186/s40364-023-00545-x
53. Bhat AA, Nisar S, Singh M. Cytokine- and chemokine-induced inflammatory colorectal tumor microenvironment: emerging avenue for targeted therapy. Cancer Commun (Lond). 2022;42(8):689–715. doi:10.1002/cac2.12295
54. Schlechter BL, Stebbing J. CCR5 and CCL5 in metastatic colorectal cancer. J IMMUNOTHER CANCER. 2024;12(5):e008722. doi:10.1136/jitc-2023-008722
55. Liu C, Yao Z, Wang J. Macrophage-derived CCL5 facilitates immune escape of colorectal cancer cells via the p65/STAT3-CSN5-PD-L1 pathway. CELL DEATH DIFFER. 2020;27(6):1765–1781. doi:10.1038/s41418-019-0460-0
56. Zhou L, Jiang Y, Liu X. Promotion of tumor-associated macrophages infiltration by elevated neddylation pathway via NF-kappaB-CCL2 signaling in lung cancer. ONCOGENE. 2019;38(29):5792–5804. doi:10.1038/s41388-019-0840-4
57. Charoentong P, Finotello F, Angelova M. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. CELL REP. 2017;18(1):248–262. doi:10.1016/j.celrep.2016.12.019
58. Fu J, Li K, Zhang W. Large-scale public data reuse to model immunotherapy response and resistance. GENOME MED. 2020;12(1):21. doi:10.1186/s13073-020-0721-z
59. Jiang P, Gu S, Pan D. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. NAT MED. 2018;24(10):1550–1558. doi:10.1038/s41591-018-0136-1