In this research, scRNA-seq data for 23 COAD tumor samples and 8 normal samples were sourced from the GSE132465 database available on the GEO website. The training and validation cohort included RNA expression data and corresponding clinical details for COAD from the TCGA database and GEO datasets GSE39582 (https://portal.gdc.cancer.gov/). Furthermore, a scrutiny encompassed 2634 ubiquitination-associated genes (URGs) sourced from the GeneCards database, each demonstrating a correlation score greater than 3.

The analysis commenced with a differential assessment to discern the variances in gene expression between tumor and normal samples in TCGA, subsequently culminating in the generation of a heatmap and a volcano plot. Subsequently, we combined the COAD samples from TCGA and GSE39582 into a merged cohort. The function normalize between arrays in R was employed to remove batch effects, and data were converted to a log2 scale before the analysis. The overall survival (OS) of the merged cohort was analyzed using a univariate Cox regression, identifying prognostically significant URGs at a p-value of less than 0.05.

Utilizing gene expression profiles, 1017 COAD samples were stratified into distinct molecular subtypes. This classification was conducted using the ConsensusClusterPlus package, which incorporates a K-means clustering algorithm to organize the samples into robust clusters, with the maximum number of clusters set at nine (maxK = 9). Employing the cumulative distribution function (CDF) curve analysis and the CDF delta area curve, the identification of the optimal number of ubiquitination-associated subtypes was undertaken, in addition to the generation of a consensus matrix heatmap to visually depict cluster affiliations. Further analysis was undertaken to explore the spatial distribution and relationship among the identified subtypes. Incorporated within this analysis were PCA and UMAP, furnishing a graphical portrayal of subtype dispersion across the samples. Additionally, differential expression of URGs and survival disparities among the subtypes were investigated using a heatmap for gene expression visualization and the ‘survival’ package for survival analysis, respectively.

Utilizing R packages “clusterProfiler” and “org.Hs.eg.db,” the analysis of differentially expressed genes across various subtypes involved leveraging Kyoto Encyclopedia of KEGG, GO enrichment, and GSEA (24). Furthermore, ssGSEA was utilized to quantify enrichment scores related to immune cell infiltration and immune functions.

Subjects with comprehensive clinical data were randomly divided into training and testing cohorts at a 1:1 ratio. Utilizing the R package “survival,” a univariate Cox regression analysis was conducted to pinpoint genes associated with prognosis (25). These genes formed the foundation for constructing a prognostic model linked to ubiquitination, employing Lasso and multivariate Cox regression methodologies (26). Post-construction, individuals were categorized into high- and low-risk groups depending on their median risk scores. Each individual’s risk score was calculated using a specific formula: R i s k   s c o r e = ∑ i = 1 n β i * e x p i , where expi represents the expression level of each URG, and βi denotes the respective gene coefficient within the signature. The Kaplan-Meier method was then used to examine overall survival (OS) differences between the two risk groups. Furthermore, the signature’s predictive accuracy was appraised with the ROC curve. Additionally, multivariate Cox regression analyses were performed in both cohorts to verify the risk score’s role as an independent prognostic marker (27).

The creation of a nomogram involved integrating the risk score, age, and pathological stage as independent prognostic variables for evaluating the likelihood of overall survival (OS) at 1, 3, and 5 years (28). To appraise the nomogram’s predictive precision, we employed the ROC, calibration, and cumulative hazard curves.

Utilizing the CIBERSORT algorithm, we computed the proportions of immune infiltrating cells in each COAD sample. Based on this data, we assessed the differences in immune cell expression across various risk groups, analyzed the correlations among immune cells, and examined their relationships with risk scores. Additionally, comparisons of tumor microenvironment (TME) scores were conducted using the R package “estimate” (29). The activation of immune checkpoints between the two groups was visualized using a barplot.

We utilized the ‘oncoPredict’ package to compute the IC50 values for various chemotherapy drugs across the two patient risk groups, aiming to gauge their sensitivity to chemotherapy. Variations between subtypes were analyzed using the Wilcoxon test.

The cell types annotated within the GSE132465 dataset were derived from prior studies. We conducted quality control on the scRNA-seq data using the “Seurat” and “SingleR” R packages. We included cells that had less than 10% of mitochondrial gene expression, more than 200 overall genes, and genes that were expressed in at least three cells, with expression levels ranging from 200 to 7000, to ensure the data's quality remained high. Subsequent analyses were conducted using the “Seurat” R package. We identified the top 2000 highly variable genes (HVGs) and employed the top 15 principal components in conjunction with these HVGs (30). Dimensionality reduction and visualization were achieved through UMAP to classify each cell type. Distinct cell types, such as T cells, B cells, and epithelial cells, were identified based on marker genes (31, 32).

We employed CellChat to assess the primary signaling inputs and outputs across all cell clusters, utilizing CellChatDB.human as the reference database. Subsequently, we applied the netVisual_circle function to illustrate the relative strength of cell-cell communication networks among the different cell clusters.

The tissue specimens were provided by Anqing First People’s Hospital affiliated with Anhui Medical University and preserved at -80°C. Ten tissue pairs, including tumor tissue (T) and precancerous tissue (N), were collected from CORD patients who underwent colon tumor resection between May 2023 and May 2024. This experiment was performed according to the previous research. All primers, along with their precise sequences, were provided by Tsingke Biotech (Beijing, China) and are detailed in Supplementary Table 1 .

For the assessment of cell proliferation, we employed the CCK-8; Vazyme, Nanjing, China). Cells were seeded at a density of 5×10³ cells per well in 96-well plates. Subsequently, the plate underwent a two-hour incubation in the absence of light at 37°C with 10 μl of CCK-8 labeling reagent per well. Cell viability was evaluated by measuring the absorbance at 450 nm using an enzymatic label reader (A33978, Thermo, USA) at intervals of 0, 24, 48, 72, 96, and 120 hours (33).

Upon reaching 95% confluency, the transfected cells were transferred to 6-well plates. A sterile pipette tip (200 μL) was used to create a straight line, followed by gentle rinsing with PBS to eliminate unattached cells and debris. The serum-free medium was then substituted to sustain cell culture. Images were taken at 0 and 48 hours in the identical position (34).

After transfecting 1000 cells, we incubated them in 6-well plates for about 14 days. At the end of this period, the cell clones became visible to the naked eye. Subsequently, the cells underwent a 15-minute fixation in 4% paraformaldehyde (PFA). Following this, staining with Crystal Violet (Solarbio, China) was conducted for 20 minutes, and the cells were air-dried at room temperature before being counted per well.

The flowchart of our study is shown in Figure 1 . We first screened 581 differential expressed genes (DEGs) from a set of ubiquitination-related genes Supplementary Figures 1 , 2 . Using 137 ubiquitination-related prognostic genes derived from univariate Cox regression analysis, consensus clustering was employed to identify the optimal number of ubiquitination-related subtypes Supplementary Figure 3 . The cumulative distribution function (CDF) diagram indicated that the CDF curve stabilized at K = 2, reflecting a robust clustering consistency. Additionally, the area under the CDF curve demonstrated notable slope alterations beyond the K values of 2 and 3. Consequently, K = 2 was selected as the most appropriate number of clusters ( Figure 2A ). Subsequently, 1017 colon cancer samples were segregated into two distinct clusters. Cluster A comprised 706 samples, while Cluster B included 311 samples. Dimensionality reduction techniques such as PCA and UAMP confirmed that the clusters occupied divergent positions and were distinctly separable, affirming the reliability of the clustering results ( Figure 2B ). Examination of gene expression and survival disparities between the clusters showed that ubiquitination-related genes were markedly upregulated in Cluster A ( Figure 2C ). Survival analysis further highlighted a significant difference in survival outcomes between the two subtypes (P < 0.001), with Cluster B exhibiting poorer survival ( Figure 2D ).

Gene Ontology (GO) enrichment analysis was conducted on the DEGs, revealing their involvement in several biological processes (BP) such as nuclear division and organelle fission. Regarding cellular components (CC), they were primarily associated with the collagen-containing extracellular matrix and chromosomal region. Concerning molecular function (MF), the differentially expressed genes (DEGs) were predominantly linked to extracellular matrix structural constituent and glycosaminoglycan binding ( Figures 3A, B ). Importantly, the activity of glycosaminoglycan binding has been reported to be associated with the prognosis of colon adenocarcinoma. Concerning KEGG analysis, we found that differential genes are mainly enriched in cell cycle and focal adhesion pathways ( Figures 3C, D ). Additionally, Gene Set Enrichment Analysis (GSEA) was conducted on both clusters, leading to the identification of enriched pathways in cluster A, such as cell cycle, DNA replication, and spliceosome ( Figure 3E ). A comparison of hallmark pathway gene signatures between the two clusters revealed distinct patterns. In contrast, cluster B exhibited enrichment in ECM-receptor interaction, focal adhesion, and cell adhesion molecule pathways as the top three signatures ( Figure 3F ). The single sample gene set enrichment analysis (ssGSEA) scores of immune cells in group A and group B were compared. The ssGSEA scores in group B were higher, and the immune infiltration and risk scores were consistent ( Figure 3G ).

Initially, a univariate Cox proportional hazards regression analysis was conducted on a merged dataset from TCGA and GEO, identifying 137 ubiquitination marker genes significantly associated with overall survival (OS, P<0.05). These genes formed the basis of a prognostic signature. The most beneficial prognostic genes were selected using a LASSO Cox regression model, and a final set of 12 genes was established through multivariate Cox regression analysis ( Figures 4A, B ). A heatmap of model gene expression is shown in Supplementary Figure 4 . The risk score was then calculated using the formula: “ (-0.815)*USP26 + (-0.193)*MYC + 0.345*OGT + (-0.408)*PRMT1 + 0.546*SNAI1 + 0.566*RPS17 + (-0.45)*RPN2 + (-0.394)*ACACA + 0.419*RNF112 + 0.503*ASNS + 0.333*MC1R + 0.46*FBXO39”. This scoring algorithm was also applied to the combined validation and overall datasets. Patients were categorized into high- and low-risk groups based on the median risk score. Notably, cluster B, as previously defined, also exhibits a higher risk score in this model ( Figure 4C ). Kaplan-Meier survival curves illustrated that individuals in the high-risk category experienced significantly worse overall survival (OS) compared to those in the low-risk group (P<0.05) ( Figures 4D–F ). The areas under the ROC curves for 1, 2, and 3 years were 0.697, 0.722, and 0.732 for the training cohort; 0.702, 0.679, and 0.673 for the testing cohort; and 0.698, 0.698, and 0.700 across all cohorts, respectively ( Figures 4G–I ). The independence of the risk score from clinical characteristics was confirmed through multivariate Cox regression analyses, as depicted in Figure 4J . The Multi-Cox regression analysis unveiled that solely age and the risk score stood as autonomous variables, wherein the hazard ratio for the risk score was determined to be 1.180 (P<0.05). The concordance index curve also underscored that the risk signature provided better predictive accuracy than other clinical features ( Figure 4K ). Furthermore, a three-dimensional scatter plot of PCA analysis illustrated the ability of the 12-gene Prognostic Risk Score (PRS) to discriminate COAD samples effectively, indicating its superior discriminatory power ( Figure 4L ).

In our effort to elucidate distinctions between risk groups, we analyzed clinical data from our merged cohort. A chi-squared test revealed significant differences in the distribution of COAD cohorts across age groups (P < 0.05) between varying risk groups, alongside notable variations in TNM stages ( Figure 5A ). Older age and more advanced TNM stages typified the high-risk cohort. To illustrate, 17% of individuals within the high-risk category were designated as stage M1, in contrast to only 10% in the low-risk subset ( Figures 5B–E ).

The construction of a nomogram aimed to assess patient risk by amalgamating clinical particulars with risk classification. Figure 5F provides a comprehensive overview of patient characteristics, including gender, age, T, N, and M stages, and risk classification. This predictive tool enables a more accurate estimation of patient risk and aids in formulating tailored therapeutic strategies ( Figure 5G ). The cumulative hazard curve indicates that patients with elevated nomogram risk exhibit a greater hazard ( Figure 5H ). To rigorously assess the predictive accuracy of the nomogram, a prognostic ROC analysis was performed, demonstrating superior performance compared to other clinical models and risk scores. The AUC values were 0.802, 0.789, and 0.760 at 1, 3, and 5 years, respectively ( Figures 5I–L ).

The figure illustrates that we used the CIBERSORT method to quantify and compare immune cell infiltration levels between the two risk groups employing the Wilcoxon test ( Figure 6A ). The violin plot presented in the study clearly demonstrates that there is a significant increase in the proportion of resting CD4 memory T cells within the low-risk group when compared to their counterparts in the high-risk group. This finding suggests a potential immune profile that is more favorable in the low-risk cohort, highlighting the importance of resting CD4 memory T cells in this context. Conversely, the high-risk group exhibited markedly higher levels of both M2 macrophages and activated mast cells, with statistical significance indicated by a p-value of less than 0.05 ( Figure 6B ). These observations can provide valuable insights into the differing immune landscapes present in varying risk categories. Furthermore, Figure 6C elaborates on the relationship between the risk score and the relative abundance of several immune cell types. It was observed that the risk score exhibited an inverse correlation with the presence of activated CD4 memory T cells, resting CD4 memory T cells, and regulatory T cells, indicating that as the risk score increases, the abundance of these beneficial T cell populations tends to decrease. In contrast, the risk score showed a positive correlation with the levels of eosinophils, neutrophils, and M2 macrophages, suggesting that higher risk scores are associated with an increase in these cell types. This differential distribution of immune cells based on risk scores points toward underlying mechanisms that may influence disease progression and patient outcomes.. In the model, the expression of the ASNS gene is inversely correlated with naive B cells, potentially indicating that the ASNS gene modulates changes in TME by regulating B cell proliferation ( Figure 6D ). We explored the relationship between risk scores and frequently identified immunotherapy biomarkers in the combined cohort. The analysis revealed that nearly all immune checkpoint genes (ICGs), including CD28 and CD70, were significantly overexpressed in the high-risk group ( Figure 6E ). The ESTIMATE methodology was employed to evaluate immune infiltration among various risk cohorts. The corresponding Figure 6F substantiated the preceding study, demonstrating that the high-risk cohort exhibited elevated estimate, stromal, and immune scores in comparison to another group.

After processing and refining the data, gene expression profiles from 53844 cells of 12 COAD samples were gathered for further analysis. Employing dimensionality reduction and log-normalization, 26 distinct cell clusters were identified ( Figure 7A ). The characterization of cell types within each cluster was achieved by comparing differentially expressed genes with canonical markers ( Figure 7B ). Differential expression of marker genes was employed to differentiate between various cellular groupings, as illustrated in Figure 7C . We then studied the ASNS gene in the model in detail. In Figures 7D, E , it is evident that the ASNS gene is highly expressed in the epithelial cells of the tumor sample. In Figure 7F , epithelial cells and stromal cells have strong cellular communication, suggesting that ASNS genes may be involved in the interaction between epithelial and stromal cells to regulate the development of colon adenocarcinoma.

Additionally, by utilizing the GDSC database, we forecasted the responsiveness of 198 drugs concerning the two risk cohorts. Within this analysis, 52 drugs exhibited varying degrees of sensitivity based on the risk stratification of the cohorts ( Figures 8A–L ). Drugs that are more sensitive to high-risk groups and have therapeutic significance are selected and shown in Figure 8 . Studies have found that potent and selective CDK9 inhibitors target transcriptional regulation in triple-negative breast cancer (35). Mitoxantrone and gemcitabine are effective in the treatment of metastatic breast cancer (36). JQ-1 (carboxylic acid), BET bromine domain inhibitors have a strong killing effect on triple-negative breast cancer cells (37). Topotecan is useful in lung cancer (38).

We knock down the ASNS gene in both SW620 and RKO cells. Figure 9A shows the transfection efficiency of the two cell types. To verify the expression level of ASNS in 10 pairs of tissues, the results showed that ASNS was highly expressed in colon cancer tissues (Figure 9B). Following in vitro testing, we gained additional insights into the function of ASNS. The CCK-8 study observed a notable decrease in proliferative activity in ASNS knockdown cells (Figures 9C, D). Similarly, the healing and migration ability of the examined cell lines were notably diminished following ASNS knockdown (Figures 9E, F). Moreover, the colony formation experiments indicated a significant reduction in the proliferation capacity of colon cancer cells after ASNS knockdown (Figures 9G, H).