Datasets were downloaded from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) using the keywords “Colorectal Cancer” or “Ulcerative Colitis.” Our data analysis process is demonstrated in Figure 1. Detailed information for each dataset, including microarray platform, sample groups, accession numbers, and sample sizes—was recorded. Only datasets containing colon tissue samples from patients with colorectal cancer and ulcerative colitis were included. A total of 6 datasets, namely, GSE52060 (Medico et al., 2013), GSE87211 (Hu et al., 2018), GSE90627 (Guo et al., 2017), GSE36807 (Montero-Meléndez et al., 2013), GSE53306 (Zhao et al., 2015), and GSE13367 (Bjerrum et al., 2010), were integrated for this study. The training set were selected as GSE52060, GSE87211, GSE36807, and GSE53306. The testing set were selected as GSE90627 and GSE13367. The details for all datasets are presented in Supplementary Table S1. To correct for batch effects from different studies, we used the “ComBat” function in the “sva” package (version 3.5.0) (Johnson et al., 2007; Leek et al., 2012). The effectiveness of batch correction was evaluated by comparing data quality before and after adjustment using principal component analysis (PCA) (Jolliffe and Cadima, 2016).

To identify key genetic alterations associated with UC and CRC, we performed differential gene expression analysis between case and control groups using the Linear Models for Microarray Data (limma) package in R. Limma is a widely used statistical tool that applies linear models to gene expression data while leveraging empirical Bayes methods to moderate the standard errors of estimated log-fold changes. This approach enhances the stability of statistical inference, particularly in studies with small sample sizes (Ritchie et al., 2015). To determine significantly differentially expressed genes (DEGs), we utilized the eBayes function, which computes moderated t-statistics, F-statistics, and log-odds of differential expression for each gene. Genes were considered significantly differentially expressed if they met the threshold of a false discovery rate (FDR) < 0.05 (adjusted p-value <0.05) and demonstrated an absolute fold change (FC) greater than 0.585 (|log₂FC| > 0.5). These stringent criteria helped ensure the robustness and reliability of our findings, highlighting genes with substantial expression changes that may play critical roles in UC and CRC pathogenesis.

To explore functional gene relationships and identify disease-associated modules, we performed Weighted Gene Co-expression Network Analysis (WGCNA). This method constructs gene co-expression networks and detects modules of highly correlated genes, often linked to specific biological traits (Langfelder and Horvath, 2008). As a crucial preprocessing step to ensure a scale-free network topology, we determined the optimal soft-thresholding power (β), selecting a β value where the scale-free topology fit index (R²) exceeded 0.8. A minimum module size of 60 genes was set to identify meaningful gene clusters.

Next, the adjacency matrix was transformed into a Topological Overlap Matrix (TOM), which enhances network robustness by reducing the effects of noise and spurious correlations. To identify gene clusters, we calculated the TOM-based dissimilarity measure (1 - TOM) and applied hierarchical clustering to group genes with similar expression patterns into modules. To refine module detection, dynamic tree cutting was implemented to segment the clustering dendrogram. To assess biological relevance, we correlated module eigengenes (principal components of modules) with clinical traits of UC and CRC. Modules with the strongest correlations and lowest p-values were selected for further analysis, helping identify key gene clusters involved in disease mechanisms and potential therapeutic targets.

A comprehensive list of cellular senescence-associated genes was obtained from the CellAge database. By intersecting the gene sets from WGCNA modules, DEGs, and the CellAge dataset via “ggvenn” package (v 0.1.9) (Gao et al., 2024), we extracted a subset of genes that are not only involved in cellular senescence but also exhibit differential expression and co-expression patterns in UC and CRC. These intersecting genes were considered as potential senescence-related biomarkers and therapeutic targets for further analysis.

Gene Ontology (GO) provides a structured, dynamically updated vocabulary encompassing gene product attributes across all species, in which GO enrichment contained 3 parts: biological processes (BPs), cellular components (CCs) and molecular functions (MFs) (Zhao et al., 2022). Kyoto Encyclopedia of Genes and Genomes (KEGG) integrates genomic, chemical, and systemic functional information, offering insights into the network of molecular interactions in the cells. For the purpose of understanding candidate genes’ function as well as participating pathways, “clusterProfiler” package (v 4.7.13) was employed for GO and KEGG analysis (Wu et al., 2021). Utilizing GO and KEGG pathway analyses, we systematically explore the functional and interactive networks that characterize the senescence landscape in UC transitioning into CRC.

We developed the final predictive model with optimal performance by applying 10-fold cross-validation on the training set, evaluating 113 model combinations derived from 12 machine learning algorithms. These algorithms included Lasso, Ridge, Stepwise GLM (Stepglm), Random Forest (RF), XGBoost, Elastic Net (Enet), Linear Discriminant Analysis (LDA), Partial Least Squares Regression for Generalized Linear Models (plsRglm), Generalized Boosted Regression Models (GBM), Naive Bayes, GLMBoost, and Support Vector Machine (SVM). The 113 models consisted of 22 individual algorithms and 91 combined algorithms, as detailed in Supplementary Table S2. To determine the best-performing model, we calculated the concordance index (C-index) for each model and selected the one with the highest C-index as the optimal model. The genes identified by this model were considered candidate disease-related genes, potentially serving as biomarkers for UC and CRC.

After constructing the integrated machine learning (IML) model, we assessed its classification performance using confusion matrices for the training set and two independent validation datasets, GSE13367 and GSE90627. To further validate the model’s predictive capability, we generated Receiver Operating Characteristic (ROC) curves for both the training and validation sets and computed the Area Under the Curve (AUC) with 95% Confidence Intervals (CI). A model was deemed statistically rational only if the ROC AUC exceeded 0.7 for both the training and validation sets (Qin et al., 2023; Chen B. et al., 2024). This approach ensured the robustness and generalizability of our model in distinguishing disease-associated genes and validating their diagnostic potential.

Differential gene expression analysis was performed using experimental data from the GEO datasets. To compare the expression levels of disease-related genes between the UC and CRC validation cohorts, we conducted Student’s t-test. Genes exhibiting statistically significant differential expression (p < 0.05) were identified as cellular senescence-related genes in UC or CRC. To evaluate their diagnostic potential, we generated Receiver Operating Characteristic (ROC) curves for each gene and calculated the Area Under the Curve (AUC) with 95% Confidence Intervals (CI). Genes with an AUC greater than 0.7 in both UC and CRC patients were considered to have significant diagnostic value (Liu et al., 2024). Furthermore, the significantly differentially expressed genes were integrated into a combined diagnostic model, and its ROC curve was constructed. If the combined model exhibited an AUC higher than that of any individual gene, it was considered a more effective diagnostic tool. The volcano plot was redrawn to visualize the upregulation or downregulation of genes with significant expression differences between UC and CRC.

A protein-protein interaction (PPI) network was constructed to explore the functional relationships and interaction dynamics among the genes with significant expression differences identified in IML. GeneMANIA (http://genemania.org/) incorporates data from multiple interaction types, including co-expression, physical interactions, genetic interactions, co-localization, pathway participation, and shared protein domains, providing a holistic view of the gene interactions. The genes with significant expression differences were input into GeneMANIA to generate a comprehensive PPI network.

To investigate disparities in immune infiltration between the two risk groups, the infiltration abundance of 22 distinct immune cell types (Chen et al., 2019) was first quantified using the CIBERSORT algorithm as implemented in the IOBR package (v 0.99.9) (Zhang et al., 2024). A Wilcoxon test was then applied to identify immune cell populations displaying significant differences (p < 0.05) between the risk groups. Subsequently, Spearman correlation analyses were performed with the psych package (v 2.4.3), using thresholds of |cor| > 0.3 and p < 0.05, to elucidate the correlation network among these differentially abundant immune cells. In addition, correlations between these immune cells and prognostic genes were evaluated under the same thresholds to further characterize the interplay between immune infiltration and gene expression profiles.

We conducted a comprehensive analysis of the expression patterns of five candidate genes (ABCB1, CXCL1, TACC3, TGFBI, and VDR) in colorectal adenocarcinoma (COAD) using publicly available RNA sequencing data from The Cancer Genome Atlas (TCGA) database. Microsatellite instability (MSI), a marker indicative of genomic instability and immunogenicity in colorectal cancer, was evaluated by computing Pearson correlation coefficients between gene expression and MSI scores. These correlations, along with significance levels (p-values), were visualized using radar charts, providing intuitive insights into the association between gene expression and MSI status (Lin et al., 2020). To examine potential differences in gene expression across tumor stages (Stages I, II, III, and IV), expression levels were visualized through boxplots, with statistical significance between stage groups determined via Wilcoxon rank-sum tests, a non-parametric method suitable for small or unevenly distributed clinical cohorts (Liu et al., 2023).

To explore potential therapeutic agents targeting cellular senescence-related genes in UC and CRC, we conducted a comprehensive drug enrichment analysis using Enrichr (https://maayanlab.cloud/Enrichr/). Initially, significant candidate compounds were screened with strict criteria, specifically applying thresholds of p-value <0.05 to ensure statistical robustness. Subsequently, we performed molecular docking analysis to validate and refine the candidate selection using the CB-Dock2 platform, an advanced version of the CB-Dock server optimized for protein-ligand blind docking. CB-Dock2 integrates cavity detection, molecular docking, and homologous template fitting to provide precise predictions of binding sites and affinities between proteins and ligands (https://cadd.labshare.cn/cb-dock2/index.php).

The selection of final candidate drugs followed a clearly defined, stepwise filtering process. Initial drug enrichment analysis identified compounds significantly interacting with the target genes. Candidate drugs passing enrichment thresholds with the lowest p-value underwent molecular docking analyses. Docking scores (Vina scores) obtained from CB-Dock2 were employed, with lower scores indicating stronger binding affinity and better potential therapeutic efficacy. Compounds with the lowest Vina scores were prioritized as potential therapeutic candidates based on their binding strength and interaction specificity, then we used CB-Dock to visualize the docking result.

Through this combined approach of enrichment analysis and molecular docking, we systematically and rigorously identified promising candidate drugs, thereby enhancing the potential for targeted therapeutic strategies in UC and CRC.

All diseased samples from the training set (GSE52060, GSE87211, GSE36807, and GSE53306) were merged into “Treat”, and all healthy control samples were merged into “Control”. In the heatmap, the validation group and experimental group are divided into different modules, and the samples in each dataset are segmented into different squares (Figure 2A). The colors of the squares represent the changes in gene expression, with red representing upregulation and blue representing downregulation. A total of 3,446 DEGs were identified by comparing the BD and control groups. Among all DEGs, 1716 genes displayed upregulation, whereas 1730 genes were downregulated (Figure 2B).

Next, WGCNA was used to identify the significant module genes associated with UC and CRC. We selected the optimal soft-thresholding power (β) to establish a scale-free topology network, ensuring that the scale-free topology fit index (R²) exceeded 0.8 (Figure 2C). The chosen β value was set to maintain the network’s scale-free characteristics. The grey module and brown module did not successfully cluster the genes commonly considered irrelevant or uninformative (i.e., the “junk module”). The turquiose (r = 0.5, p = 8 × 10⁻³⁰) module displayed the highest correlation with UC and CRC (Figure 2D). The relationship between module membership and gene significance in the turquiose module is calculated (Cor = 0.93, p < 10^–200) and plotted (Figure 2E). Consequently, 3,571 significant module genes were identified.

Through a comprehensive analysis integrating data from the CellAge database, weighted gene co-expression network analysis (WGCNA) modules, and differentially expressed genes (DEGs), we identified a refined subset of 112 shared genes implicated in cellular senescence and their association with UC and CRC. The Venn diagram (Figure 2F) summarizes the intersection results, illustrating the overlap between the gene sets and emphasizing the genes that could serve as pivotal links between cellular senescence and disease progression in UC and CRC.

The differentially expressed genes identified using the limma R package were analyzed through Gene Ontology (GO) enrichment analysis (Figures 3A,B). The results were ranked in ascending order based on adjusted p-values (p.adjust) and GeneRatio. In the Biological Processes (BP) category, the top three pathways with the lowest p.adjust values and the highest number of enriched genes were morphogenesis of a branching structure, morphogenesis of a branching epithelium, and gland development. In the cellular components (CC) category, the top 3 terms were cytoplasmic vesicle lumen, secretory granule lumen, and collagen−containing extracellular matrix. In the molecular functions (MF) category, the top 3 terms were DNA−binding transcription factor binding, DNA−binding transcription activator activity, DNA−binding transcription activator activity, RNA polymerase II−specific. The KEGG pathway analysis revealed key pathways that were significantly enriched among the genes identified in our study (Figures 3C,D). These pathways included the PI3K-Akt signaling pathway, p53 signaling pathway, and the cell cycle, which are known to play pivotal roles in regulating cellular senescence, survival, proliferation, and apoptosis. The involvement of these pathways underscores the potential mechanisms through which cellular senescence could influence the transition from UC to CRC.

A comprehensive machine learning approach involving 12 algorithms was implemented with a 10-fold cross-validation process to identify the most robust diagnostic model based on 112 shared genes (Figure 4A). The analysis was conducted using the training dataset and validated across two external datasets (GSE90672 and GSE13367). The final model, which demonstrated the best performance, was constructed by integrating Stepglm [both] and Enet [α = 0.6]. Specifically, the Stepglm [both] algorithm identified 10 pivotal genes, including ABCB1, AGR2, BCL2L1, CXCL1, FOXO1, SOX4, TACC3, TGFβI, VDR, and VEGFA, while the Enet [α = 0.6] algorithm optimized the model’s reliability. The validation datasets remained completely independent and were not involved at any stage of feature selection, model training, parameter tuning, or optimization, thereby preventing any potential data leakage or information contamination. Furthermore, all cross-validation procedures, feature selection steps, and modeling approaches were performed exclusively within the training set. The calibration curves, illustrated in Figures 6D,E, show high AUC values for the training set (Figure 4B, AUC = 0.991), as well as the testing set GSE90627 (Figure 4C, AUC = 1.000) and GSE13367 (Figure 4D, AUC = 0.993), indicating a strong agreement between the predicted probabilities and observed clinical outcomes. These results highlight the robust calibration and diagnostic performance of the proposed model.

Ten pivotal genes were included in the following ROC analysis. All 10 pivotal genes showed high significance (p < 0.001) in CRC and its control group (Figure 5A), while ABCB1, CXCL1, TACC3, TGFβI, VDR displayed high significance (p < 0.001) in UC and its control group (Figure 5B). Based on the significant differences in gene expression, ABCB1, CXCL1, TACC3, TGFβI, and VDR were integrated into a combined model. All the 5 genes were included in ROC analysis. We calculated the AUC values for each gene and the combined model separately (Figure 5C). The results showed that the AUC values of all genes were not less than 0.7, and the AUC value of the combined model (AUC = 0.989) was higher than that of any individual gene. Therefore, the combined model has greater diagnostic value compared to any individual gene.

The PPI network for genes which were included in combined model genes (ABCB1, CXCL1, TACC3, TGFβI, and VDR) were created via the GeneMANIA database (https://genemania.org/). In the GeneMANIA map, a total of 20 genes (CXCL5, CXCL6, POSTN,SLC22A3, SLC22A1, CCL11, TACC2, TACC1, CKAP5, APCS, ACKR1, ACBC4, RXRB, MEDI, CYP3A4, SLC22A2, NDEL1, CSCR2, BAG1, and CLIP4) were found to have gene interactions with five combined model genes (Figure 6). In the GeneMANIA network, physical interactions between pivotal genes and other genes account for 77.64%, while co-expression accounts for 8.01%, demonstrating the strong protein-protein interactions within the GeneMANIA network topology.

Immune correlation analysis was performed with all samples in training set (Figure 7A). The infiltration landscape showed that 22 kinds of immune cell distributions in the control and treat groups. Fourteen types of immune cells (neutrophils, mast cells activated, mast cells resting, macrophages M2, macrophages M0, monocytes, NK cells activated, T cells follicular helper, T cells CD4 memory activated, T cells CD4 memory resting, T cells CD4 naive, T cells CD8, and B cells memory) infiltrated significantly (p < 0.001) between the control and treat groups (Figure 7B). Correlation analysis between immune cells indicates that Macrophage M2 exhibited significantly negative correlation with activated T cells CD4 naive (r = −0.63, p < 0.05), T cells CD8 had positive correlation with macrophages M2 (r = 0.31, p < 0.05) (Figure 7C).

The correlation between genes and 22 immune cell types, as well as the interrelationships among immune cells, has been systematically analyzed and visualized (Figure 7D). Among the findings, ABCB1 exhibits the strongest positive correlation with T cells CD4 memory resting, while demonstrating the most pronounced negative correlation with Monocytes (Figure 7E). Similarly, CXCL1 is most positively correlated with Neutrophils and Macrophages M0, whereas its most significant negative associations are observed with NK cells activated and Macrophages M2 (Figure 7F). In the case of TACC3, its highest positive correlation is identified with Neutrophils, whereas its strongest negative correlations are noted with B cells memory and Plasma cells (Figure 7G). Notably, TGFβI shows the most significant positive correlation with Mast cells activated, while displaying a marked negative correlation with Plasma cells, T cells CD8, T cells gamma delta, and Macrophages M2 (Figure 7H). The gene VDR did not exhibit strong correlations with immune cells in the analysis (Figure 7I).

ABCB1 expression significantly decreased with advancing tumor stage, with higher median expression observed in early-stage (Stage I/II) compared to late-stage tumors (Stage III/IV), consistent with earlier reports of its downregulation in colorectal carcinogenesis (Figure 8A). Conversely, CXCL1 and TACC3 were progressively upregulated in advanced stages. CXCL1 showed significantly elevated expression in Stage III–IV tumors compared to Stage I–II, indicative of enhanced inflammation (Figure 8B). Similarly, TACC3 levels significantly increased in late-stage disease, aligning with its previously reported association with tumor progression and poor prognosis in colorectal cancer (Figure 8C). TGFBI and VDR exhibited no clear stage-dependent expression patterns. TGFBI expression fluctuated without significant differences, while VDR levels remained relatively stable across stages, aligning with previous studies reporting limited stage-dependent variation (Figures 8D,E). MSI analysis revealed significant positive correlations for TACC3 and CXCL1 (p < 0.05), suggesting their potential roles in MSI-high tumor biology, whereas ABCB1, TGFBI, and VDR lacked significant associations with MSI status (Figure 8F).

Enrichr database was utilized to screen therapeutic agents targeting the five combined model genes. The analysis highlighted several compounds potentially effective in targeting genes associated with UC and CRC, and the compounds (p < 0.05) are listed in Supplementary Table S2. Protein sequences were obtained from Uniprot (https://www.uniprot.org/), and compound structures were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/). The identified therapeutic agents included iodoquinol, cefaclor, pyrithione, 5-Aminosalicylic acid, 2-Mercaptobenzothiazole, 1,10-Phenanthroline, eugenol, Bisulfite, Alitretinoin, and gossypol (Figures 9A,B).

Genes CXCL1 (Uniprot ID: P09341), TACC3 (Uniprot ID: Q9Y6A5), and VDR (Uniprot ID: P11473) were selected as target proteins, and their interactions with candidate drugs were visualized through interconnected curves, enabling clearer exploration of drug-target relationships (Figure 9C). The four most significant drugs, iodoquinol (CID: 3,728), cefaclor (CID: 51,039), pyrithione (CID: 1570), and 5-Aminosalicylic acid (CID: 4075), were chosen based on the lowest enrichment p-values for further validation using molecular docking. Docking results, represented by Vina scores, are summarized in Table 1. Lower Vina scores indicate stronger binding affinities, justifying the prioritization of these compounds as therapeutic candidates. Visualization of docking interactions was performed using PyMOL software (Figures 9D–K). Importantly, 5-Aminosalicylic acid is already clinically established for UC treatment, providing additional validation and supporting the overall reliability and validity of our model.