Overview of the expression and genetic variation of ARGs in CRC

To investigate the impact of ARGs on the occurrence and progression of CRC, we analyzed the RNAseq data of normal and tumor tissues from the TCGA database. We identified DEGs using a threshold of FDR < 0.05 and an absolute logFC greater than 1. This analysis led to the identification of 50 differentially expressed ARGs, which were visualized through a heatmap and a volcano plot (Figure S2A, B). We further performed univariate Cox analysis on the DEGs to determine their association with overall survival in CRC patients. The results revealed that a total of 13 differentially expressed ARGs were significantly associated with patient survival outcomes, with 8 ARGs associated with unfavorable survival outcomes (Figure S2C). The prognostic relevance and interrelationships of these ARGs were summarized in a prognostic-related network depicted in Figure S2D. Additionally, we analyzed copy number variations (CNVs) to further explore ARG mutations in colorectal cancer. We determined the occurrence rate of ARG mutations by integrating CNV data and identified 45 ARGs with CNV alterations (Figure S2E), mostly characterized by copy number amplifications. The genomic locations of these CNV alterations are shown in Figure S2F. These results shed light on the potential prognostic value of ARGs in colorectal cancer and demonstrate their association with patient survival outcomes. The analysis of CNV alterations provides additional insights into the occurrence and distribution of mutations in ARGs, which may contribute to a deeper understanding of the molecular mechanisms underlying CRC development.

Identification of ARGs subtypes

We employed unsupervised cluster analysis to stratify CRC patients’ samples based on ARGs expression levels, leading to the identification of two distinct clusters, termed Cluster A and Cluster B (Figure S3). Survival analysis indicated a favorable survival advantage for patients in Cluster B (Fig. 1A). Subsequently, we applied three common dimensionality reduction techniques, namely PCA, t-SNE, and UMAP, to visualize the data, with each cell population represented by different colors. The visualizations clearly distinguished Cluster A and Cluster B (Fig. 1B-D). Differential expression analysis identified 13 ARGs with significant expression differences between the two clusters, including 9 upregulated and 4 downregulated in Cluster A compared to Cluster B (Fig. 1E). Notably, the clinical characteristics also exhibited significant differences between the two subtypes, with patients in Cluster A displaying older age and higher STAGE staging (Fig. 1F).

Fig. 1

Molecular subtypes based on ARGs in CRC and their clinicopathological features. (A) Kaplan-Meier survival analyses for the two molecular subtypes. UAMP(B), tSEN(C), and PCA(D) presented a great difference between the A and B subtypes. (E)The expression levels of prognosis related differentially expressed ARGs between A and B subtypes. (F) The heatmap showed the prognosis related differentially expressed ARGs expression profiles and clinicopathologic characteristics among subtypes A and B

To explore the immune microenvironment of Clusters A and B in CRC, we performed ssGSEA to compare the expression profiles of 23 immune cell subtypes between the two clusters. The results revealed that 18 out of 23 immune cell subtypes exhibited differential expression between the clusters (Fig. 2A). Furthermore, to gain insights into the underlying biological alterations associated with the distinct clusters, we conducted GSEA and GSVA analyses. GSEA unveiled significant activation of ECM-receptor interaction and Focal adhesion pathways in Cluster A, whereas ribosome and oxidative phosphorylation-related metabolic pathways were enriched in Cluster B (Fig. 2B). GSVA analysis of the two subtypes indicated that highly expressed DEGs in Cluster A were significantly enriched in cell communication and information pathways, MAPK signaling pathway, and cell adhesion-related pathways (Fig. 2C, D).

Fig. 2

Tumor microenvironment and Functional enrichment analysis of ARGs-based clusters in CRC. (A) Analysis of infiltrating immune cells between A and B subtypes. (B) Analysis of GSEA for subtypes A and B. C, D. Analysis of GSVA for subtypes A and B

Metabolic reprogramming induced by ARGs subtypes

The metabolic pathway covers various important biological molecule synthesis and decomposition processes within cells, which are crucial for normal cell function and survival. At the same time, dysmetabolism is also a core feature of cancer. The further screening of metabolism-related pathways through GSEA reveals that the molecular subtypes of CRC based on ARGs also reflect the metabolic alterations involved in the progression of colorectal cancer. As shown in Additional Table 1, and 2, Besides fatty acid metabolism, Cluster B exhibits significantly enriched metabolic pathways, encompassing the following categories: Carbohydrate Metabolism (Butanoate Metabolism, Propanoate Metabolism, Pyruvate Metabolism, CITRATE CYCLE TCA CYCLE, Starch and Sucrose Metabolism), Nucleotide Metabolism (Pyrimidine Metabolism and Purine Metabolism), Amino Acid Metabolism (Arginine and Proline Metabolism, Cysteine and Methionine Metabolism), Drug Metabolism (Drug Metabolism-Other Enzymes and Drug Metabolism-Cytochrome P450), Vitamin Metabolism (Retinol Metabolism), Organic Acid Metabolism (Glyoxylate and Dicarboxylate Metabolism), and Other Metabolic Pathways (Porphyrin and Chlorophyll Metabolism, Glutathione Metabolism).

Construction and validation of risk prognostic models

To establish a feature-scoring model for evaluating the role of ARGs in CRC, we employed LASSO regression analysis to select the best prognostic feature-related genes from prognostic-related key genes. After incorporating the variables into the LASSO regression model with minimized λ, five ARGs were selected to construct the risk-scoring model (Fig. 3A). Initially, the samples were stratified into train and test groups with a 1:1 ratio. Figure 3B-D depicted the gene expression levels and survival time, along with the status distribution between the high and low-risk groups in various groups. The results showed an increasing proportion of patient mortality with elevated risk scores. Simultaneously, survival analysis across the three groups revealed that patients with higher risk scores exhibited significantly inferior overall survival (OS) compared to those with lower risk scores (Fig. 4A-C). Time-dependent ROC analysis for 1-year (train group 0.695, test group 0.677, all samples group 0.850), 3-year (train group 0.850, test group 0.814, all samples group 0.817), and 5-year (train group 0.813, test group 0.788, all samples group 0.771) overall survival further validated the robust predictive capacity of the ARGs-associated risk model for colon cancer patient survival (Fig. 4D-F). Subsequent multivariable Cox regression results indicated that the risk score independently served as a prognostic factor for OS (Fig. 4H). Patients in the A and B clusters also demonstrated differences in risk scores, with Cluster A patients exhibiting higher risk scores, corresponding to a poorer prognosis for previous patients in Cluster A (Fig. 4I). Sankey plots illustrating both groups displayed the correlation between ARGs risk model grouping, ARGs subtypes, and survival status. The high-risk group displayed a higher proportion of fatal outcomes compared to the proportion of alive patients, and within Cluster A, the proportion of high-risk patients was higher than that of low-risk patients (Fig. 4J). Based on the results of univariate and multivariate Cox regression analyses, we constructed a nomogram incorporating clinical staging and ARGs (Fig. 5A). The cumulative graph demonstrated a significant distinction between high and low-risk groups over time (Fig. 5B). Calibration curves indicated the high accuracy of the nomogram (Fig. 5C). These findings collectively suggest that ARGs possess a reliable capacity to discriminate tumor outcome differences.

Fig. 3

Construction of prognostic model based on ARGs in CRC (A) LASSO Cox regression analysis. Distribution of the heatmap of ARGs (upper), survival time (middle), and risk score (below) in all cohorts (B), training cohort(C), and test cohort(D)

Fig. 4

K-M survival curves of low- and high-risk patients in the training cohort(A), test cohort(B), and all cohort(C). D-F. The ROC curves at 1-, 3- and 5-year in the mentioned three cohorts. H. Forest plots showing the results of the multivariate Cox regression analysis. I. Differences in risk scores between clusters A and B J. Sankey diagram shows the relationship between different ARG scores, risk scores, and survival outcomes

Fig. 5

Construction and validation of a nomogram (A) Nomogram for predicting the 1-, 3-, and 5-year OS of CRC patients. (B) Draw a cumulative risk map for high and low-risk groups (C) Calibration plots show the fits of 1-, 3- and 5-year predictions

Tumor mutation burden and immune landscape of ARGs subtypes

The current investigation employed the Maftools software package to undertake a comparative analysis of somatic mutations between high and low-risk oncology cohorts. Our analysis revealed that the high-risk cohort (Figure S4A) displayed a significantly higher tumor somatic mutation rate when compared to the low-risk cohort (Figure S4B). Further identification of tumor driver genes in both cohorts indicated a higher number of such genes in the low-risk group (Figure S4C, D). The analysis of co-occurring and mutually exclusive mutations in the high and low-risk cohorts showcased the correlation among mutated genes, with green and orange colors representing co-occurrence and mutual exclusivity, respectively (Figure S4E, F). Moreover, the forest plot illustrating differentially mutated genes between the high and low-risk cohorts (Figure S5) indicated a higher proportion of mutated genes in the low-risk group.

Immunogenomic patterns and immunotherapy analysis

Using the CIBERSORT method, we determined the proportions of 22 immune infiltrating cell types in each sample and subsequently investigated the differential expression of immune cells between the high and low-risk groups. Our findings demonstrated that the high-risk group exhibited a higher proportion of resting immune cells (Fig. 6A). Furthermore, the correlation analysis of the 22 immune cell types revealed positive correlations indicated in red and negative correlations indicated in blue (Fig. 6B). Leveraging these results, we examined the relationship between genes involved in model construction and immune cells, observing that most genes showed negative correlations with immune cells (Fig. 6C). Interestingly, we observed a higher immune functionality in the high-risk group (Fig. 6D).

Fig. 6

Tumor microenvironment features related to the ARGs-based signature in CRC. (A) The proportions of immune cells in each sample of low- and high-risk groups according to CIBERSORT analysis. (B) The correlation of immune cells (C) The heatmap showed the relationship between ARGs and immune cells. (D) Immune function analysis

To evaluate the responsiveness of patients in the high and low-risk groups to immunotherapy, we initially conducted a risk score analysis and explored the correlation of immune checkpoint molecules. Remarkably, we identified ten immune checkpoint molecules, including ICOS, HAVCR2, CTCN1, etc. which displayed significant positive correlations, while five immune checkpoint molecules, including ICOSLG, TNFSF9, IDO2, CD40LG, and TNFRSF4, showed significant negative correlations (Fig. 7A, B). Additionally, ESTIMATE analysis revealed that the high-risk group exhibited higher immune, stromal, and ESTIMATE scores compared to the low-risk group (Fig. 7C). Subsequently, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) to assess the efficacy of immunotherapy in high and low-risk populations. The results depicted in Fig. 7D indicated a higher proportion of immunotherapy-resistant patients in the high-risk group. Furthermore, the high-risk group displayed lower microsatellite instability (MSI) scores, higher TIDE scores, and T cell exclusion scores, with no significant difference in T cell dysfunction compared to the low-risk group (Fig. 7E). Finally, through the evaluation of PD1 and CTLA-4, we found that patients in the low-risk group may exhibit a more favorable response to immunotherapy (Fig. 7F).

Fig. 7

Immune checkpoint profiles and immunotherapy evaluation related to the ARGs-based signature in CRC. A, B. correlations between the immune checkpoint expression and risk score. C. Comparison of TEM scores between the high and low-risk groups. D. analyze immunotherapy responses in high and low-risk groups by TIDE E. Predictive immunotherapy by comparing the score of MSI, Exclusion, Dysfunction, and TIDE in high and low-risk groups. F. Immunotherapy response analysis of high and low-risk groups through TCIA

Drug response analysis

In order to identify potentially effective drugs for treating colorectal cancer, we examined the correlation between drug sensitivity (IC50 values) and gene expression profiles in high-risk and low-risk groups. Our analysis revealed that several drugs exhibited promising effects on colorectal cancer patients. Specifically, 39 drugs were identified as having potential efficacy in the treatment of colorectal cancer. Some of these drugs include Afuresertib, AGI-5198, Crizotinib, Dabrafenib, Dasatinib, and many others. Our analysis revealed that certain drugs, such as Dasatinib, Doramapimod, JQ1, and NU7441, demonstrated greater efficacy in the low-risk subgroup of patients (Figure S6). On the other hand, the remaining drugs showed higher efficacy in the high-risk subgroup of patients. These findings suggest that different patient risk profiles may influence the response to specific drugs, highlighting the importance of personalized medicine approaches in tailoring treatment strategies for colorectal cancer patients. Further investigation is warranted to understand the underlying mechanisms and potential biomarkers associated with drug response in different risk groups.

Validation of key genes expression in tumor cell lines in vitro

We detected the expression of 5 ARGs (Cd36, FABP4, ANGPT1, ACADL, and CPT2) in CRC cell lines (SW480 and HCT116) and normal intestinal epithelial cells (NCM460) using RT-qPCR in vitro. Except for ANGPT1, the expression trends of the other ARGs were consistent with those in TCGA (Figure S7).