RNA-seq data and clinical details for COAD patients were sourced from TCGA database (https://portal.gdc.cancer.gov/). A total of 585 transcriptomic datasets, clinical data (GSE40967), and single-cell datasets from 13 COAD samples (GSE110009) were acquired from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). A total of 576 ARGs were sourced from the GeneCards database (https://www.genecards.org/) and the Harmonizome database (https://maayanlab.cloud/Harmonizome/) (19). To maintain the accuracy and reliability of the analyses, only ARGs with a correlation coefficient exceeding 0.4 were chosen for further examination.

In this research, we obtained transcriptomic and clinical data from COAD samples within TCGA database. Differential expression analysis was performed using the R package “DESeq2” to identify differentially expressed genes (DEGs) with a fold change of at least 2 and a p-value of less than 0.01. By intersecting these DEGs with ARGs, we identified 134 significantly different genes in COAD. Standardized merging with GEO data led to further prognostic analysis, revealing 37 genes associated with overall survival (OS). Gene copy number variations were assessed using data obtained from the UCSC Xena website (https://xena.ucsc.edu/), while a protein-protein interaction (PPI) network was generated using the online tool STRING to investigate co-expression relationships and potential molecular interactions among the genes (20, 21).

We performed an in-depth clustering analysis of samples using the K-means clustering algorithm from the “ConsensusClusterPlus” package, setting the maximum number of clusters to 9 (maxK = 9). Through detailed analysis of the consistency matrix of clustering results, we successfully determined the optimal number of clusters to be 2, which was validated using principal component analysis (PCA) (22). Heatmap analysis demonstrated relationships between alterations in gene expression and clinical characteristics, while survival curve analysis showed that Cluster A had a better prognosis than Cluster B, revealing notable differences in immune cell expression between the two clusters. Gene Set Variation Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA) were utilized to further investigate biological functions and pathway activity differences across the various subtypes.

We utilized the “createDataPartition” function in R to randomly partition the dataset into training and testing subsets. Univariate Cox regression analysis was performed to identify significant genes, followed by Lasso regression and cross-validation to optimize a multivariate Cox model (23). The “glmnet” package facilitated the calculation of risk scores for COAD patients based on the refined model. Patients were categorized into high-risk and low-risk groups based on the median risk score. The risk score was calculated using the formula: Risk score = ∑(expi * βi), where expi represents the expression level of each gene and βi denotes the corresponding regression coefficient. Kaplan-Meier survival curves were plotted using the “survival” package, and time-dependent receiver operating characteristic (tROC) curves were generated with the “timeROC” package to evaluate the model’s predictive accuracy regarding patient survival.

Nomograms were developed using the “rms” and “regplot” packages in R to forecast survival in COAD patients, taking into account factors such as age, sex, and stage. Calibration curves validated the precision of the nomogram’s predictions for survival rates at 1 year, 3 years, and 5 years (24). Decision curve analysis revealed that the nomogram’s predictive ability was superior to that of individual clinical factors. Performance assessment of the risk score showed that its accuracy in prediction surpassed that of conventional clinical indicators.

Both univariate and multivariate Cox regression analyses indicated a significant association between age, stage, and risk scores with the survival duration of patients with COAD. The risk score model constructed using these variables demonstrated significant differences in survival times across patients of varying gender, age, and stage.

Using R packages “clusterProfiler” and “org.Hs.eg.db,” we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on differentially expressed genes across different risk groups. GO enrichment analysis revealed distinct patterns in the distribution of differentially expressed genes between high-risk and low-risk groups, spanning biological processes (BP), molecular functions (MF), and cellular components (CC) (25). Additionally, KEGG pathway analysis highlighted the key pathways enriched by these differentially expressed genes.

We employed the CIBERSORT algorithm for precise calculation of immune cell infiltration proportions (26). A heatmap was utilized to visually represent the strength and directionality of correlations among different immune cell types. Specific immune cell subpopulations exhibited significant expression differences across risk groups. Furthermore, we employed the “estimate” package to calculate tumor microenvironment scores (TME scores) to evaluate the infiltration levels of immune and stromal cells within the tumor microenvironment (27).

We performed thorough quality control and filtering of scRNA-seq data utilizing the “Seurat” and “SingleR” R packages. Each gene was required to be expressed in a minimum of three cells and to have an expression level of at least 50 genes. The “subset” function was applied to filter cells based on the criteria of having more than 50 genes and less than 5% mitochondrial gene expression. Data normalization was carried out using the “NormalizeData” function with a scaling factor of 10,000. The “FindVariableFeatures” function was then employed to identify genes exhibiting high variability, selecting the 1,500 genes with the most significant expression fluctuations for further analysis. Through PCA dimensionality reduction and t-SNE clustering analysis, we successfully identified multiple distinct cell clusters and recognized several marker genes (28). We used the SingleR algorithm in R to annotate cell types in our single-cell dataset by comparing it to a reference dataset. The results were visualized using t-SNE, which provided a clear representation of the cellular landscape. Additionally, we applied the Monocle algorithm to infer the developmental trajectories of the cells and construct dynamic models of cellular differentiation.

We conducted an extensive drug screening utilizing the gene prediction drug online resource (https://maayanlab.cloud/Enrichr/), selecting drugs based on an adjusted p-value of less than 0.05. The 2D structures of these drugs were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) and converted into 3D models using “Chem3D.” Protein structure data for the relevant genes were sourced from the PDB database (http://www.rcsb.org/). The receptor was prepared using “PyMol,” removing water molecules and small ligand compounds. Molecular docking was conducted to identify active pockets using “Autodock Vina v.1.5.7” (29).

The cell lines present in this study were obtained from The Cell Resource Center of Shanghai Life Sciences Institute. Firstly, lipo3000 transfection reagent was mixed with siRNA to form transfection complex. These complexes are then added to cells in a petri dish to bind to the cell membrane and enter the cell. After transfection, the cells were cultured for 24 hours under suitable culture conditions so that the foreign nucleic acid could be expressed or function. Finally, PCR was used to detect the transfection effect and the knockdown of exogenous genes. Primer information and siRNA sequence are Supplementary Table 1 .

Total protein from the cells was extracted using RIPA Lysis Buffer (Beyotime, P0013B), and its concentration was determined using the Enhanced BCA Protein Assay Kit (Beyotime, P0010). The results were measured by ImageJ software.

After diluting the cells to an appropriate concentration, they are evenly seeded into sterile 6-well plates and then cultured in an incubator for 10-14 days until visible cell colonies form. Next, the cells are fixed with formaldehyde and stained with crystal violet.

Cultivate the cells to near confluence and then gently create a straight line scratch in the cell layer at the bottom of the culture dish using a sterile pipette. Subsequently, remove any cell debris from the scratch area and rinse once with fresh medium to eliminate floating cells. After that, place the cells back in the incubator and observe and capture images of the scratch healing at 0 and 48 hours to assess the rate and capacity of cell migration.

In this study, we performed logarithmic transformation and batch correction to standardize the data. All data analyses and graphical visualizations were performed using R software (version 4.3.3). To explore the co-expression relationships among genes, we employed Pearson correlation analysis. Furthermore, Spearman correlation analysis was performed to assess the association between risk scores and levels of immune cell infiltration. A comprehensive predictive model was established by combining univariate, Lasso, and multivariate Cox regression analyses. For comparative evaluations, t-tests were utilized to determine differences between two groups, while one-way ANOVA was applied to examine statistical variations across multiple groups. A p-value of less than 0.05 was deemed significant, with * denoting P < 0.05, ** indicating P < 0.01, and *** representing P < 0.001.

Our research workflow is clearly illustrated ( Figure 1 ). We collected transcriptomic data and clinical details from 483 COAD samples and 41 adjacent normal samples via TCGA database. After data organization, we merged the samples and converted probe IDs into gene IDs. We then compared gene expression levels between tumor and normal tissues, performing differential analysis using the “DESeq2” package (fold change = 2, p-value = 0.01). The intersection with ARGs yielded 134 significantly different genes in COAD ( Figure 2A ). We then standardized and merged the transcriptomic dataset of 585 samples and corresponding clinical data from the GEO database with TCGA dataset, extracting the expression levels of DEGs. Employing survival status and OS as outcome measures, univariate Cox regression analysis identified 37 prognostic genes, which are depicted in a forest plot ( Figure 2B ). The prognostic network diagram depicted in Supplementary Figure 1A illustrates the co-expression relationships among 29 high-risk and 8 low-risk genes. We subsequently displayed the frequency of copy number increases or decreases for specific genes within the samples, indicating the prevalence of gene variations across the samples, as shown in Supplementary Figure 1B . A copy number circle plot depicted in Supplementary Figure 1C shows the frequency of gene copy numbers at corresponding chromosomal positions. We constructed a PPI network for the anoikis-related DEGs, as shown in Supplementary Figure 1D , identifying genes with protein interaction relationships through connecting lines. The correlation network diagram in Supplementary Figure 1E illustrated the strength of correlations among different genes.

Based on the expression levels of ARGs, the samples were classified into two subtypes, A and B ( Figure 2C ). PCA demonstrated that the subtyping effectively separated the samples into distinct groups ( Figure 2D ). To further investigate the differences between subtypes A and B, survival curves confirmed a significant disparity in patient survival (p < 0.001), with subtype A exhibiting superior prognosis compared to subtype B ( Figure 2E ). The subtype heatmap integrated subtype classifications with clinical characteristics, providing a detailed analysis of the distribution of upregulated and downregulated genes ( Figure 2F ). We also examined the differential expression of ARGs and immune cells between subtypes A and B ( Figures 2G, H ). Subsequently, GSVA and GSEA analyses were performed to investigate the biological functions and pathways associated with the different subtypes, emphasizing the discrepancies in biological pathway activities ( Figure 3A ). The findings indicated that cell adhesion molecules (CAMs), cytokine-cytokine receptor interactions, ECM receptor interactions, focal adhesion, and neuroactive ligand-receptor interactions were significantly reduced in subtype A ( Figure 3B ). Conversely, subtype B exhibited increased expression of CAMs, cytokine-cytokine receptor interactions, ECM receptor interactions, neuroactive ligand-receptor interactions, and systemic lupus erythematosus ( Figure 3C ).

We integrated a transcriptome dataset of 1,034 samples with survival data from TCGA and GEO databases, randomly partitioning it into a training cohort (n=517) and a test cohort (n=517). We employed the Least Absolute Shrinkage and Selection Operator (LASSO) to construct a regression model, identifying the point of minimal error through cross-validation ( Figures 4A, B ). Subsequently, we performed multivariate Cox regression analysis to enhance the model, selecting nine genes as risk features: NAT1, INHBB, FGF2, CD36, CCDC80, SPP1, MMP3, S100A11, and GZMB ( Figure 4C ). The risk score for COAD patients was computed using the following formula: Risk Score (RS)=(−0.25374×NAT1 expression)+(0.19775×INHBB expression)+(0.14497×FGF2 expression)+(0.23801×CD36 expression)+(−0.34139×CCDC80 expression)+(0.09980×SPP1 expression)+(−0.10291×MMP3 expression)+(0.35502×S100A11 expression)+(−0.16830×GZMB expression) Using the median risk score, we stratified the samples into high-risk and low-risk groups within both the training and test cohorts. We evaluated our model’s predictive efficacy for OS in COAD patients across all cohorts. Kaplan-Meier survival analysis revealed that high-risk COAD patients experienced significantly poorer OS in all three cohorts ( Figures 4D–F ). The ROC curves for the three cohorts showed AUC values exceeding 0.6, specifically 0.687, 0.709, and 0.705 for 1-year, 3-year, and 5-year survival, respectively ( Figures 4G–I ). These results indicate that our model is highly accurate in predicting survival outcomes for COAD patients at 1, 3, and 5 years.

We constructed a nomogram ( Figure 5A ) that aggregates scores from each clinical characteristic (age, sex, stage) to predict patient survival. To assess the predictive performance of the nomogram, we generated calibration curves and decision curves. The calibration curves showed high accuracy in forecasting 1-year, 3-year, and 5-year survival using the nomogram ( Figure 5B ). The decision curves for 1 year, 3 years, and 5 years indicated that the nomogram’s predictive power significantly exceeded that of other clinical factors ( Figures 5C–E ). In analyzing the survival rates of COAD patients, we compared the effectiveness of ARG features. Our findings revealed that the risk score exhibited the highest C-index, affirming its predictive accuracy and surpassing traditional clinical indicators, including pathological stage, age, and sex ( Figure 5F ).

Univariate Cox regression analysis was conducted to generate a forest plot ( Figure 6A ), showing significant associations between survival time and age, stage, and risk score. Multivariate Cox regression analysis further demonstrated statistical relationships between survival time and age, T, M, N stages, and risk score, as illustrated in a forest plot ( Figure 6B ). Consequently, the risk score calculated based on the nine ARGs can predict the survival rate of COAD patients. Stratification by age, sex, and stage demonstrated significant differences in survival times between high-risk and low-risk groups across all categories (p < 0.01). High-risk patients showed reduced survival durations compared to those classified as low-risk ( Figures 6C–H ).

To examine the distribution patterns of gene sets across BP, MF, and CC, we conducted GO enrichment analysis using differentially expressed genes from the high and low-risk groups. GO analysis indicated significant enrichment in cellular response to chemokine, endopeptidase activity, and immunoglobulin complex ( Figures 7A ). KEGG pathway analysis revealed that the differentially expressed genes were mainly enriched in pathways including cytokine-cytokine receptor interaction, IL-17 signaling pathway, viral protein interactions with cytokines and receptors, chemokine signaling pathway, and rheumatoid arthritis ( Figures 7B ). Notable differences in expression levels were found for plasma cells, activated CD4 memory T cells, resting NK cells, monocytes, M1 macrophages, and M2 macrophages across the different risk categories ( Figure 7C ). We conducted a differential analysis of immune cells between high-risk and low-risk groups, illustrated by a heatmap that shows the correlation strength and direction among various immune cell types, such as monocytes, eosinophils, and activated dendritic cells ( Figure 7D ). A heatmap illustrating the correlations between immune cells, model genes, and risk scores was also presented ( Figure 7E ). Notably, risk scores exhibited significant expression differences with M0 macrophages, M1 macrophages, M2 macrophages, resting NK cells, plasma cells, activated CD4 memory T cells, and CD8 T cells, with the strongest correlation identified with activated CD4 memory T cells (R = -0.4, p < 2.2e−16) ( Figures 7F–L ). Patients in the low-risk group generally exhibit higher TME scores, which may be associated with stronger immune responses and lower tumor purity, thereby potentially having a positive impact on patient treatment responses and survival outcomes. The violin plot demonstrated statistically significant differences in tumor TME scores, including immunescore and estimatescore, between high and low-risk groups ( Figure 7M ).

We retrieved single-cell datasets for COAD from the GEO database, utilizing gene expression profiles from 13 COAD samples for further analysis. Initial quality control and filtering were performed, isolating the 1,500 most variably expressed genes, with the top ten most variable genes highlighted ( Figure 8A ). PCA was employed to reduce dimensionality on these 1,500 genes ( Figure 8B ), and expression levels of PC1-PC4 feature genes were visualized in a heatmap ( Figure 8C ). TSNE clustering analysis of 15 principal components identified 17 distinct clusters ( Figure 8D ). The differential analysis for each cluster has identified 5,933 marker genes, which are detailed in Supplementary Table 2 . We conducted a detailed examination of the expression patterns for the nine key genes within the constructed model. Notably, the INHBB gene exhibited significantly elevated expression in cluster 15, while CCDC80 showed higher expression in cluster 13, and S100A11 was notably expressed in cluster 12 ( Figures 8E, F ). Analysis of cell types within each cluster revealed diverse populations, including T cells, B cells, epithelial cells, monocytes, macrophages, fibroblasts, tissue stem cells, and natural killer cells ( Figure 8G ). Further annotation revealed specific gene-cell type associations: INHBB was most highly expressed in endothelial cells, CCDC80 in fibroblasts, and S100A11 in epithelial cells. These findings provide crucial molecular evidence for understanding the roles of different cell types in biological processes. The cells were primarily derived from hematopoietic stem cells and differentiated along two distinct pathways: one pathway leading to the formation of epithelial-like cells, and the other leading to the generation of monocytes ( Figure 8H ).

In the prognostic model, the INHBB gene was identified as the most significant biomarker affecting OS in COAD patients, as shown in Supplementary Figure 2 . To explore potential therapeutic strategies targeting this key gene, we conducted a comprehensive drug screening using a gene-drug prediction database. By considering adjusted p-values, we successfully identified 22 compounds with potential therapeutic relevance in Supplementary Table 3 . To further assess the binding affinity of these candidate drugs with their targets, we performed molecular docking for the top four compounds. Using AutoDock Vina v.1.5.7, we analyzed the binding sites and interactions of these candidates with the INHBB protein, calculating binding energies for each interaction. The docking results indicated that INHBB had the lowest binding energy with risperidone at -9.5 kcal/mol, suggesting a stable interaction. Visualize the docking results of four drugs with INHBB using Pymol software ( Figures 9A–D ). This analysis not only allows us to predict potentially effective drugs but also provides critical insights into drug-target interactions, informing drug development and optimization efforts.

Further analysis revealed that patients with high INHBB expression in COAD had worse OS and disease-free survival (DFS) ( Figures 10A, B ). In RKO and SW620 cell lines, INHBB knockdown resulted in significant changes in INHBB expression levels ( Figure 10C ). The verification of the low efficiency protein level is shown in the Supplementary Figure 3 . In the wound healing assay, we observed a significant decrease in the proliferation activity of RKO and SW620 cells post-INHBB knockdown compared to control cells ( Figures 10D, E ). To further validate the impact of INHBB on the proliferative capacity of colorectal cancer cells, we conducted colony formation assays. The results showed that after INHBB gene knockdown, both cell lines exhibited reduced colony number and size ( Figure 10F, G ).