Single-cell transcriptomics reveals the prognostic and immune infiltration significance of FOLFOX-bevacizumab treatment and lysine crotonylation characteristics in colon cancer

Introduction

Colorectal cancer (CRC) ranks third in worldwide malignancy incidence and second in mortality. By 2040, the CRC burden is projected to increase to 3.2 million new cases and 1.6 million deaths (1,2). Colon Cancer (CC), a malignancy arising from the colonic mucosal epithelium and an aggressive CRC subtype, accounts for approximately two-thirds of CRC cases compared to rectal cancer (2,3). Its development is closely linked to genetic variations, environmental factors, and lifestyle. Despite prominent advances in CC screening, early detection and prevention remain challenged by tumor biological heterogeneity and dynamic progression (4). Traditional CC treatments include surgery, chemotherapy, and radiotherapy; nevertheless, outcomes vary considerably and are often associated with severe side effects.

FOLFOX (oxaliplatin, leucovorin, and 5-fluorouracil)-bevacizumab (FOLFOX-Bev) is a targeted therapeutic regimen for CC. Studies indicate that in advanced CC, FOLFOX-Bev demonstrates optimal efficacy against cancer cells exhibiting high proliferation, metastasis, and pro-angiogenic properties (5). However, resistance is prevalent (6). Meanwhile, FOLFOX-Bev has varying therapeutic effects on patients with different genotypes. The adenomatous polyposis coli (APC) mutation status of patients affects bevacizumab treatment (7), and its cellular resistance mechanism needs further in-depth analysis.

Recent research highlights the central role of epigenetic modifications in tumor progression. Lysine crotonylation (Kcr), was initially identified as a reversible protein modification on lysine residues within histones enriched at promoters and enhancers in human cells and male germinal cells (8). It involves the transfer of a four-carbon unsaturated acyl group (CH₃CH=CHC=O⁻) to the ε-amino group of target protein lysines, serving crotonyl-CoA as the donor, catalyzed by crotonyl transferases (e.g., p300/CBP, MOF) (9). Kcr is implicated in various diseases and biological processes, including cancer, tissue injury, and inflammation (10). Some studies report elevated Kcr modification levels in CC tissues compared to adjacent normal tissues (11). However, how Kcr modification influences FOLFOX-Bev treatment response in CC at the single-cell level has not been systematically elucidated.

Overall, this study integrated single-cell and The Cancer Genome Atlas (TCGA) datasets to analyze the dynamic changes in the CC microenvironment and the functional role of Kcr modification pre-/post-FOLFOX-Bev treatment, exploring their associations with immune infiltration, cell communication, and malignant phenotypic transformation. Subsequently, a prognostic model was constructed and validated employing least absolute shrinkage and selection operator (LASSO)/multivariate Cox regression to assess differential immunotherapy responses, predict suitable drugs for different risk groups, and identify relevant pathways involving differentially expressed genes (DEGs) as potential therapeutic targets, thereby offering novel strategies for precision therapy. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-703/rc).

Methods

Data collection

The dataset of 452 CC patients was downloaded from TCGA (https://portal.gdc.cancer.gov/) as the training set.

Utilizing the GEO (https://www.ncbi.nlm.nih.gov/), the single-cell dataset (GSE232525) of one CC patient pre-/post-FOLFOX-Bev treatment and the data of 55 CC patients (GSE17537), as the validation dataset, were obtained.

Genes related to “Kcr” were searched and downloaded from the GeneCards Human Gene Database (https://www.genecards.org/) (Relevance score >2.0), yielding 102 genes (Table S1).

Gene expression signature matrices for tumor-infiltrating immune cell (TIIC) were downloaded from the GSVA platform (http://software.broadinstitute.org/gsea/msigdb/index.jsp).

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Analysis and annotation of single-cell RNA sequencing (scRNA-seq) data

scRNA-seq data (GSE232525) were preprocessed via the “Seurat” R package (12). Step 1: cells were retained if mitochondrial gene percentage was <50%, unique genes detected in a single cell (nFeature_RNA) were between 500 and 10,000, and total RNA molecules detected in a single cell (nCount_RNA) were between 200 and 6,000. Step 2: data were normalized by implementing “NormalizeData”, “FindVariableFeatures”, and “ScaleData” functions in the R package to identify the top 2,000 highly variable genes. Step 3: batch correction and dimensionality reduction: batch effects between samples were removed using the function “FindIntegrationAnchors” in the R package with CCA. Based on the elbow plot (Figure S1), the function “FindNeighbors” (dims =35) was utilized to construct a SNN graph. Step 4: cells were clustered via “FindClusters” and confirmed using the “clustree” R package (resolution =0.5) (Figure S2), resulting in 21 clusters visualized through Uniform Manifold Approximation and Projection (UMAP).

Marker genes for each cluster were identified by the function “FindAllMarkers” (logfc.threshold =0.585, min.pct =0.25, only.pos=T). Based on marker genes in the literature (5) and the CellMarker database (http://xteam.xbio.top/CellMarker/), combined with the automated annotation results from the “singleR” package (13), the cells in each cluster were annotated to determine and annotate the different cell clusters. Bubble plots, violin plots for the expression levels of common marker genes, and heatmaps of the top 10 DEGs were generated.

Copy number variation (CNV) feature identification and epithelial cells subpopulation analysis

Epithelial cells were obtained through the function “subset” and reprocessed (steps 1–3 above). Subsequently, the data was dimensionally reduced by applying RunUMAP, and FindClusters was employed to identify the clusters. Eight epithelial subclusters were yielded via the “Clustree” (resolution =0.2). The inferCNV tool was implemented to infer CNV, thereby identifying malignant cells. The results were visualized to confirm malignant clusters. GSVA analysis of the hallmark epithelial-mesenchymal transition (EMT) pathway (https://www.gsea-msigdb.org/gsea/msigdb/) was performed to compare discrepancies of scores between identified malignant/non-malignant epithelial subpopulations.

Single-cell scoring and cell communication

The “AUCell” R package (14) was applied to score individual cells based on the Kcr-related gene set. Cells were classified as high/low Kcr (dividing standard =0.05) based on the threshold calculation results. The proportion of high/low cells within each major cell type was compared according to this standard.

Malignant cells were stratified into high-/low-Kcr groups based on the threshold from the calculation of area under the curve cell-level enrichment analysis (AUCell). Cell communication was analyzed between all subpopulations (including these two malignant groups) using the “CellCall” R package (15) to compare discrepancies in the two groups. Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, a ligand-receptor-transcription factor (LR-TF) dataset was collected. According to the prior knowledge of interactions, CellCall inferred intercellular communication by combining the expression of certain ligand-receptor (L-R) pairs’ ligands/receptors and the activity of downstream transcription factor (TFs).

Weighted gene co-expression network analysis (WGCNA)

The score of the Kcr-related gene sets for the TCGA training set samples was calculated through single-sample Gene Set Enrichment Analysis (ssGSEA) via the “GSVA” package (16). The function “goodSamplesGenes” in the “WGCNA” R package (17) was applied to determine whether the genes of the samples need to be filtered. Utilizing the “pickSoftThreshold” tool, the soft threshold of 7, the minimum number of module genes to 100, and the threshold to 0.2 were formed to merge the modules within the distance threshold. The module genes with the absolute value of correlation being the largest and P<0.05 were selected.

Identification of sensitive and insensitive subpopulations

Subpopulation analysis on malignant cells was conducted in the same way as for epithelial cells (resolution =0.1). Ro/e (ratio of observed to expected cell proportion) analysis was applied to identify the differences in distinct cell clusters pre-/post-treatment. Changes in the malignant subclusters pre- vs. post-FOLFOX-Bev treatment were analyzed to identify sensitive and drug-resistant subpopulations. Kcr scores were compared between the two groups. DEGs between sensitive and high-Kcr and drug-resistant and low-Kcr cells were identified [|avg_log₂ fold change (FC)| >0.25 & P_val_adj <0.05].

Prognostic model construction

The intersection of WGCNA module genes and DEGs between “sensitive & high-Kcr” and “drug-resistant & low-Kcr” groups was identified. Within the TCGA training set, the intersecting genes were further screened using LASSO regression analysis, where the optimal parameter (λ) was determined through ten-fold cross-validation. Subsequently, the selected genes were incorporated into a multivariate Cox regression analysis to construct the final prognostic model.

The risk score for each patient was calculated based on the regression coefficients (Coefficient) and gene expression values (expressionvalue) by the following formula:

Riskscore=∑Coefficient(gene)×Expressionvalue(gene)[1]

Patients in the training/validation sets were stratified into high-/low-risk groups based on the median risk score. Kaplan-Meier (K-M) survival curves were then generated using the “survival” R package (18) to compare the survival differences between the high-/low-risk groups (P<0.05). Furthermore, time-dependent receiver operating characteristic (ROC) curves were plotted to assess the prognostic performance of the model.

Prognostic nomogram construction and validation

Univariate/multivariate Cox regression analyses were performed to assess the independent prognostic value of the risk score. A prognostic nomogram integrating the risk score and clinical factors was constructed. The predictive ability of the model was evaluated through the decision curve analysis (DCA) curve and the calibration curve.

Analysis of immune microenvironment in high-/low-risk groups

The abundance of TIICs in the two risk groups was calculated by the ssGSEA algorithm from the “GSVA” R package. The TIIC gene expression signature matrix was downloaded from the GSVA platform. Moreover, the expression levels of human leukocyte antigen (HLA)-related genes and common immune checkpoint genes (19) were compared between risk groups (Wilcoxon test, P<0.05).

Prediction of immunotherapy response of patients

Tumor Immune Dysfunction and Exclusion (TIDE) (20) was applied to predict immunotherapy response in distinct subpopulations by simulating two primary mechanisms of tumor immune escape: T cell dysfunction induced by high cytotoxic T lymphocyte (CTL) levels and T cell infiltration suppression at low CTL levels.

Identification and analysis of DEGs in high-/low-risk groups

DEGs between the two risk group samples in the training dataset were identified utilizing the “edgeR” package (21) [|logFC| >1 & false discovery rate (FDR) <0.05]. Enrichment analysis (qvalue <0.05) of DEGs was performed employing the “clusterProfiler” R package (22) and “enrichplot” R package (23) to explore potential biological processes.

Tumor mutation and tumor heterogeneity analysis

Somatic mutation [single nucleotide variants (SNVs)] data from TCGA were analyzed via the “maftools” package (24) to plot mutation profiles in high-/low-risk groups. Differences in mutation types, SNV class, and mutation rates between groups were analyzed. The top 20 mutated genes were identified.

Drug prediction

Drug responses in high- and low-risk groups of the colon adenocarcinoma (COAD) dataset were predicted using the “oncoPredict” package (25). Specifically, a predictive model was built using drug sensitivity data from the Genomics of Drug Sensitivity in Cancer 2 (GDSC2) database and gene expression profiles of cancer cell lines. Tumor gene expression profiles of patients in the two risk groups of TCGA-COAD were input. We assumed that cells with similar expression profiles exhibited similar sensitivity to drugs. By comparing the patient’s expression profile with those of hundreds of cancer cell lines in GDSC2, the predicted half-maximum inhibitory concentration (IC₅₀) for each drug was predicted for each patient. Subsequently, the predicted IC₅₀ values were compared between the two risk groups to screen for drugs with significant differences (intergroup differences were defined using Wilcox.test, ns: P>0.05; *: P≤0.05; **: P≤0.01; ***: P≤0.001; ****: P≤0.0001). A correlation diagram was drawn between drug prediction IC₅₀ and gene expression levels.

Statistical analysis

All statistical analyses were performed using R software (version 4.4.2). Kaplan-Meier analysis was used to generate survival curves, and differences between groups were compared using the log-rank test. Univariate and multivariate Cox regression analyses were conducted to evaluate the independent prognostic value of the risk score. LASSO regression analysis was performed with ten-fold cross-validation to determine the optimal λ parameter. Time-dependent ROC curves were used to assess the predictive performance of the prognostic model. Differences in immune-related gene expression and drug sensitivity between the high- and low-risk groups were compared using the Wilcoxon test. Differentially expressed genes were identified using the edgeR package with the criteria of |logFC| >1 and FDR <0.05. Enrichment analyses were considered statistically significant at q value <0.05. Unless otherwise specified, P<0.05 was considered statistically significant.

Results

scRNA-seq data processing and analysis

Following data screening and analysis using the “Seurat” package, 21 cell clusters were identified and visualized via UMAP (Figure 1A). A heatmap of the top 5 most significantly DEGs per cluster was generated (Figure 1B, P<0.05). Cell type annotation identified 7 major cell types: B, natural killer (NK), T, myeloid, epithelial, endothelial cells, and fibroblasts (Figure 1C). Dot plots (Figure 1D) and violin plots (Figure 1E) were generated based on key marker gene expression.

Figure 1 Cell type identification from single-cell sequencing data. (A) UMAP dimensionality reduction plot of single-cell transcriptome, identifying 21 cell clusters. (B) Heatmap of the top 5 DEGs across all cell clusters. (C) UMAP plot of cell type annotation results. (D) Marker gene expression patterns for 7 cell types. (E) Violin plot of expression levels of canonical marker genes for each cell type. DEGs, differentially expressed genes; NK, natural killer; UMAP, Uniform Manifold Approximation and Projection.

To deeply analyze the cell heterogeneity within single-cell samples, the following analyses were conducted. First, a heatmap was drawn for the top 10 genes with the most marked differential expression among the seven main cell types (Figure 2A, P<0.05). Next, the distribution of different cell types was analyzed (Figure 2B). Finally, the proportion composition of cell types in each sample was presented (Figure 2C). The results indicated that both the untreated group (JCML1) and the drug treatment group (JCML2) contained various immune cells (myeloid, T, B, NK cells) and stromal cells (epithelial, endothelial cells, fibroblasts). The difference was that post-FOLFOX-Bev treatment, B, T, and NK cells were enriched, and the numbers of other types of cells (fibroblasts, endothelial, epithelial cells) decreased. At the same time, the reduction in myeloid cells was the most severe.

Figure 2 Cellular composition analysis of single-cell samples. (A) Heatmap of the top 10 most significantly DEGs for each cell type. (B) Quantitative distribution of different cell types. (C) Bar plot of proportions of cell types in each sample (JCML1: pre-FOLFOX-bevacizumab therapy; JCML2: post-FOLFOX-bevacizumab therapy). DEGs, differentially expressed genes; NK, natural killer.

Epithelial subpopulation analysis and malignant cell identification

To identify malignant subpopulations within epithelial cells, 8 clusters were yielded by dimensionality reduction clustering (Figure 3A). Employing myeloid cells as a reference, inferCNV analysis of CNV revealed large-scale chromosomal amplifications and deletions in epithelial cells (Figure 3B). Comparison of CNV scores demonstrated appreciably higher genomic instability in epithelial cells versus myeloid cells (P<0.0001) (Figure 3C). Analysis of CNV scores per epithelial cluster indicated Cluster 1/4 had lower scores than other groups; thus, they were defined as normal epithelial cells, while others were cancer cell clusters (Figure 3D).

Figure 3 Identification of malignant cells within epithelial cells. (A) UMAP plot of epithelial cell subpopulations (8 clusters). (B) inferCNV heatmap of genomic copy number variation (red: chromosomal amplification; blue: deletion). (C) Boxplot of CNV scores between epithelial cells and myeloid cells. (D) Violin plot of CNV scores for epithelial subclusters. (E) Classification of epithelial cell subpopulations (normal vs. malignant). (F) Distribution of malignant cells across all cell types. (G) Violin plot comparing EMT pathway scores between normal and malignant epithelial cells (****, P<0.0001). CNV, copy number variation; EMT, epithelial-mesenchymal transition; inferCNV, inference of copy number variations; UMAP, Uniform Manifold Approximation and Projection.

To verify malignant cell identification reliability, GSVA analysis based on the EMT pathway was performed. Epithelial cells were divided into two functional clusters: normal epithelial cells and malignant cells (Figure 3E). Mapping malignant cells to the whole-cell atlas presented their distribution (Figure 3F). Comparing EMT activity between malignant and normal epithelial cells revealed notably higher EMT scores in malignant cells (P<0.0001) (Figure 3G), confirming the reliability of malignant cell identification.

AUCell scoring of single-cell data

To resolve cellular heterogeneity in Kcr modification, single cells were scored through the AUCell algorithm and classified into high/low Kcr scoring groups (threshold =0.05). The distribution of high/low Kcr scoring cells aligned with the score distribution (Figure 4A,4B). Within malignant cells, the proportion of cells with high Kcr scores was lower than that with low Kcr scores, indicating that Kcr is associated with the progression of CC malignant cells. (Figure 4C). Further comparison of Kcr modification levels within the same cell types pre-/post-treatment indicated pronounced increases in malignant, epithelial, B, and myeloid cells post-treatment (P<0.05) (Figure 4D,4E). This suggests FOLFOX-Bev treatment remarkably enhances Kcr gene set activity.

Figure 4 Cellular characterization of Kcr modification. (A) UMAP plot of distribution of Kcr gene set scores. (B) Spatial distribution map of high/low Kcr scoring groups. (C) Bar plot of proportions of high/low Kcr score cells across cell types. (D) Violin plot of AUCell scores pre-/post-therapy across cell types. (E) Comparison of score changes across cell types pre-/post-therapy. AUC, area under the curve; AUCell, area under the curve cell-level enrichment analysis; Kcr, lysine crotonylation; NK, natural killer; UMAP, Uniform Manifold Approximation and Projection.

Malignant subpopulation and CellCall cell communication analyses

To resolve malignant cell heterogeneity, subpopulation clustering was performed, and Ro/e analysis revealed distribution changes pre-/post-treatment (Figure 5A). Results demonstrated that post-treatment, the proportion of malignant cells decreased, while normal epithelial cells increased (Figure 5B). In the high Kcr-scoring group, the Ro/e index increased post-treatment, indicating enrichment, while the low-scoring group showed the opposite trend (Figure 5C). As Cluster 0 presented the most significant post-treatment Ro/e decrease, it was defined as the sensitive subpopulation (sensitive), with others classified as drug-resistant (insensitive). Ro/e value changes per cluster are detailed in Figure 5D. To further validate the association between Kcr status and treatment response, we compared the proportions of high- and low-Kcr cells in sensitive and resistant malignant cells. The proportion of low-Kcr cells was significantly higher in the resistant subpopulation than in the sensitive subpopulation, while high-Kcr cells were relatively enriched in the sensitive subpopulation (Figure S3, P<0.05). This result, from a reverse grouping perspective, further supports the pronounced association between Kcr-related molecular status and FOLFOX-Bev treatment response. Next, the spatial distribution of sensitive and insensitive subpopulations pre-/post-treatment was compared. Pre-treatment, sensitive cells formed large, predominant clusters, while resistant cells were scattered as small clusters interwoven with sensitive cells. Post-treatment, sensitive cells decreased sharply and became dispersed, indicating treatment efficacy. Resistant cells increased in proportion and showed altered centralized distribution, highlighting their spatial position (Figure 5E).

Figure 5 Heterogeneity analysis of malignant cell treatment response. (A) UMAP plot of spatial distribution of malignant cells. (B) Proportional changes in cell types pre-/post-therapy. (C) Comparison of Ro/e index in high/low Kcr score groups pre-/post-therapy. (D) Changes in Ro/e values for malignant cell clusters pre-/post-therapy (+++: Ro/e >1; ++: Ro/e =0.8–1; +: Ro/e =0.2–0.8). (E) Evolution of spatial distribution of sensitive and drug-resistant cell groups pre-/post-therapy. Kcr, lysine crotonylation; NK, natural killer; UMAP, Uniform Manifold Approximation and Projection.

The CellCall tool was applied to explore intercellular communication within the microenvironment. Strong communication between T cells, fibroblasts, and endothelial cells (Figure 6A). Comparing L-R interaction strength between different cell types revealed that highly active L-R pairs in the high Kcr-scoring group included CXCL16-CXCR6, EDN3-EDNRA, and BMP2-BMPR2 (Figure 6B). Subsequently, downstream signaling pathways in cells from high/low Kcr-scoring groups were analyzed. Sankey diagrams of LR-TF axes for high/low Kcr-scoring group cells sending signaling to T cells, fibroblasts, and endothelial cells were generated (Figure 6C-6H).

Figure 6 CellCall cell communication analysis (a toolkit for cell-cell communication analysis). (A) Circular plot of intensity of sending and receiving cell communication for each cell type (malignant cells stratified into two groups by Kcr score). (B) Heatmap of L-R interaction intensity between cell types within high/low Kcr groups. (C) Sankey diagram of LR-TF axes from high-Kcr score cells to T cells. (D) Sankey diagram of LR-TF axes from low-Kcr score cells to T cells. (E) Sankey diagram of LR-TF axes from high-Kcr score cells to fibroblasts. (F) Sankey diagram of LR-TF axes from low-Kcr score cells to fibroblasts. (G) Sankey diagram of LR-TF axes from high-Kcr score cells to endothelial cells. (H) Sankey diagram of LR-TF axes from low-Kcr score cells to endothelial cells. Kcr, lysine crotonylation; L-R, ligand receptor; LR-TF, ligand receptor-transcription factor; NK, natural killer.

Prognostic model construction and validation

To further mine diagnostic genes, WGCNA was performed in the training set. With R²>0.9 and β=7 as optimal parameters (Figure 7A), the minimum module size was set to 100, the merge threshold to 0.2, and modules merged under the distance threshold, identifying 15 modules (Figure 7B). The MEturquoise (turquoise) module (Cor =0.57, P<0.05) and MEmidnightblue (midnight blue) module (Cor =0.49, P<0.05) were notably positively correlated with Kcr levels. Further calculation of module gene correlation with module eigengenes (MEs) is shown in Figure 7C,7D. The turquoise module (Cor =0.75, P<0.001) and midnight blue module (Cor =0.78, P<0.001) were markedly positively correlated with the Kcr phenotype.

Figure 7 Key gene screening via WGCNA. (A) WGCNA analysis of scale-free topological fit index (R²) and soft threshold β=7. (B) Association between Kcr and MEs. Columns represent Kcr features; rows represent consensus modules. Numbers indicate correlations between traits and module eigengenes, with P values in parentheses. (C) Correlation between GS and MM for turquoise module genes. (D) Correlation between GS and MM for midnightblue module genes. (E) Venn of intersection of WGCNA module genes and DEGs. Cor, correlation coefficient; DEG, differentially expressed gene; GS, gene significance; Kcr, lysine crotonylation; ME, module eigengene; MM, module membership; RGs, related genes; WGCNA, weighted gene co-expression network analysis.

Previous experimental results showed that after treatment, the number of malignant cells reduced, accompanied by higher Kcr levels. Therefore, we believed that sensitive and high-Kcr cells are the most easily killed cell clusters. Conversely, drug-resistant and low-Kcr cells should be the more resistant cell type. Therefore, we focused specifically on the differences between the “sensitive & highly Kcr” cell subset and the “drug-resistant & low-Kcr” cell subset. Through the comparison, 387 DEGs were identified (Table S2). Finally, to screen core genes with both Kcr regulatory and therapeutic response value in malignant cells, the WGCNA-selected MEturquoise and MEmidnightblue module genes were intersected with the 387 DEGs, yielding 70 key genes (Table S3). These genes exhibited consistent signals in multi-scale analyses, demonstrating correlations at the levels of overall Kcr phenotype, transcriptional co-expression network, and single-cell therapy response. These signals provide a foundation for constructing biologically based prognostic models (Figure 7E).

Subsequently, LASSO Cox regression analysis (λ=0.02615987) screened 5 genes from the 70 candidates: USP53, UACA, C7orf50, NOP53, IDH2 (Figure 8A,8B). Multivariate Cox regression then constructed the final prognostic model, selecting 4 key genes: USP53, UACA, C7orf50, and IDH2. USP53 and IDH2 were protective factors [hazard ratio (HR <1)], while UACA increased risk (HR >1), and C7orf50 exhibited a potential risk trend but insufficient evidence (Figure 8C). The risk score for each patient was calculated by the following formula:

Figure 8 Prognostic model construction and evaluation. (A) Coefficient distribution plot of log(λ) sequence in LASSO model. (B) LASSO Cox regression coefficient path plot. (C) Forest plot of multivariate Cox regression results. (D) Risk scores of each patient in the TCGA colon cancer set. (E) Survival outcome statistics of each patient in the TCGA set. (F) Risk scores of each patient in validation set (GSE17537). (G) Survival outcome statistics for each patient in validation set (GSE17537). (H) ROC curve of TCGA cohort. (I) Kaplan-Meier survival curves of high/low-risk groups in TCGA cohort. (J) ROC curve of GSE17537 cohort. (K) Kaplan-Meier survival curves of high/low-risk groups in GSE17537 cohort. *, P<0.05; ***, P<0.001. AIC, Akaike information criterion; AUC, area under the curve; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Riskscore = −0.82747 × USP53 + 0.79311 × UACA + 0.40645 × C7orf50 − 0.34001 × IDH2

The risk score was calculated for each patient. CC samples in the training set (TCGA) and validation set (GSE17537) were stratified into high-/low-risk groups based on the median risk score. Survival analysis (including survival status statistics and risk score distribution) indicated lower survival rates in the high-risk group (Figure 8D-8G). To assess the prognostic model’s value, K-M survival curves and ROC curves were plotted, and AUC values were calculated. In the TCGA training set, high-risk patients had lower survival rates, and the model showed good predictive performance at 1/3/5 years in the ROC analysis (AUC at 1 year =0.72, AUC at 3 year =0.67, AUC at 5 year =0.68) (Figure 8H,8I). Similar results were obtained in the validation set (GSE17537): high-risk patients had lower survival rates, and the model predictive performance at 1/3/5 years was relatively high in the ROC analysis (AUC at 1 year =0.66, AUC at 3 year =0.68, AUC at 5 year =0.7) (Figure 8J,8K).

Furthermore, both univariate and multivariate Cox regression analyses confirmed that the risk score had independent prognostic value (Figure 9A,9B). A nomogram constructed based on the risk score and clinical pathological parameters predicted the 1/3/5-year survival rates of CC patients (Figure 9C). The calibration curve and DCA results further verified the reliability of this nomogram in predicting survival rates (Figure 9D,9E).

Figure 9 Validation of prognostic model. (A) Univariate Cox regression analysis of clinical features and risk score. (B) Multivariate Cox regression analysis of clinical features and risk score. (C) Prognostic nomogram of the integration between clinical features and risk score. (D) Calibration curves of 1-, 3-, and 5-year predictions of the nomogram. (E) Decision curve analysis for 1-, 3-, and 5-year predictions of the nomogram. f-1/3/5: net benefit of the prognostic model for predicting 1-, 3-, and 5-year survival, respectively; All-1/3/5: net benefit when all patients are assumed to receive intervention at 1-, 3-, and 5-year time points, respectively; None: no intervention strategy. CI, confidence interval; DCA, decision curve analysis; M, metastasis; N, lymph node; T, tumor.

Genetic mutation analysis in high-/low-risk groups

Mutation types and frequencies were analyzed to compare genetic mutations between high-/low-risk score groups. We also analyzed and plotted mutation profiles in the TCGA dataset high-/low-risk groups through maftools. Results demonstrated both groups predominantly featured Missense_Mutation, with other mutation types being less prevalent, indicating gene mutations are common in CC and unaffected by risk score (Figure 10A,10B). Genetic changes differed between the two risk groups. In the TCGA data, the total mutation frequency of model genes was higher in the low-risk group. Specifically, multi-hit and missense mutations in UACA and frameshift and missense mutations in USP35 were relatively prominent in the high-risk group. The low-risk group more commonly featured various frameshift mutations in UACA and missense and nonsense mutations in USP35 (Figure 10C,10D). Additionally, although high-frequency mutation genes like APC, TP53, KRAS, TTN, and PIK3CA were found in both groups, their mutation rates differed: PIK3CA mutation was more common in the high-risk group (32%), while MUC16 mutation was more common in the low-risk group (30%) (Figure 10E,10F). These results revealed vital discrepancies in key driver gene mutation characteristics between the two risk groups.

Figure 10 Tumor mutation analysis. (A) Mutation profile summary of high-risk group. (B) Mutation profile summary of low-risk group. (C) Mutation status of model genes in high-risk group. (D) Mutation status of model genes in low-risk group. (E) The top 20 genes with highest mutation rates in high-risk group. (F) The top 20 genes with highest mutation rates in low-risk group. DEL, deletion; INS, insertion; ONP, oligonucleotide polymorphism; SNP, single-nucleotide polymorphism; SNV, single nucleotide variant; TMB, tumor mutational burden.

TIIC analysis

Subsequently, TIIC abundance in high-/low-risk groups was calculated utilizing the ssGSEA algorithm. Significantly increased infiltration of plasmacytoid dendritic cells (pDCs) in the high-risk group. Concurrently, neutrophil and NK cells infiltration decreased markedly. Increased pDCs infiltration was accompanied by elevated HLA and major histocompatibility complex (MHC) molecule levels, suggesting relatively stronger antigen-presenting capacity in the high-risk group. Conversely, decreased neutrophil and NK cell infiltration indicated weakened innate immune response function, potentially playing a more negative role in disease risk control within the high-risk group (Figure 11A). Additionally, immune checkpoint analysis revealed that CD40LG was lowly expressed in the high-risk group, while programmed cell death protein 1 (PDCD1) was highly expressed (P<0.05) (Figure 11B). These differences suggested CD40LG and PDCD1 may be potential immunotherapy targets. Finally, analysis of HLA-related gene abundance within cells revealed appreciably higher levels of HLA-E, HLA-F, HLA-DQB2, HLA-DRB5, HLA-A, HLA-DRB1, HLA-DPB1, HLA-C, and HLA-B in the high-risk group (P<0.05) (Figure 11C).

Figure 11 Immune characterization of high-/low-risk groups. (A) ssGSEA analysis of immune cells and functions in TCGA cohort. (B) Expression of immune checkpoint genes in high-/low-risk groups. (C) Differential expression of HLA genes between high-/low-risk groups (P values from Wilcoxon test; ns: P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001). aDCs, activated dendritic cells; APC, antigen-presenting cells; CCR, C-C chemokine receptor; DCs, dendritic cells; HLA, human leukocyte antigen; iDCs, immature dendritic cells; IFN, interferon; MHC, major histocompatibility complex; NK, natural killer; ns, not significant; pDCs, plasmacytoid dendritic cells; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; Tfh, T follicular helper cells; Th, T helper; TIL, tumor-infiltrating lymphocytes.

Immunotherapy evaluation and drug sensitivity prediction

Immunotherapy assessment revealed markedly higher TIDE scores and a lower predicted proportion of responders to immunotherapy in high-risk CC patients (P<0.05) (Figure 12A,12B), suggesting this group is less likely to benefit from immunotherapy. Using the “oncoPredict” package, 6 sensitive drugs were screened for the two risk groups (Figure 12C). For patients in the high-risk group, six drugs (TAF_5496_1732, OF_1853, Erlotinib_1168, AZD3759_1915, KRAS. G12C Inhibitor.12_1855, and Sinularin_1838) showed higher sensitivity, suggesting that these six predicted drugs may have better treatment effects in the high-risk population. Meanwhile, Dasatinib_1079, Sepantronium.bromide_1941, AZD7762_1022, AZD1332_1463, IGF1R_3801_1738, and WEHI.1997 may have better treatment effects in the low-risk group (Figure 12C).

Figure 12 Therapeutic response prediction of high-/low-risk groups. (A) Comparison of TIDE scores between high-/low-risk groups (Wilcoxon test). (B) Comparison of predicted immunotherapy response proportions (Fisher’s exact test). (C) Drug sensitivity analysis of high-/low-risk groups. ***, P<0.001; ****, P<0.0001. TIDE, Tumor Immune Dysfunction and Exclusion.

DEGs and functional annotation in high-/low-risk groups

To explore functional differences between high-/low-risk groups, differential expression analysis identified 889 DEGs (Figure 13A). Gene Ontology (GO) analysis revealed gene enrichment in cell fate regulation (e.g., Wnt signaling pathway), tissue development/barrier function (e.g., epidermis development), neural signaling (e.g., acetylcholine receptor signaling pathway), and mucosal immunity (e.g., mucosal immune response) (Figure 13B). KEGG analysis further indicated these genes are active in signal transduction (e.g., olfactory transduction), neuroactive L-R interactions, and hormone signaling pathways (e.g., hormone signaling) (Figure 13C). These results collectively point to systematic differences in cell behavior decisions, immune microenvironment, and signal transduction networks between the two groups.

Figure 13 DEGs and functional annotation of high-/low-risk groups. (A) Volcano plot of DEGs between high-/low-risk groups. (B) Bar plot of GO enrichment analysis for DEGs between high-/low-risk groups. (C) Bar plot of KEGG pathway enrichment analysis for DEGs between high-/low-risk groups. BP, biological process; CC, cellular component; DEGs, differentially expressed genes; FC, fold change; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Discussion

CC development is a complex multi-step process mediated by multiple factors, including genetic and epigenetic dysregulation under microenvironmental influences (26,27). Epigenetic modifications are known to regulate pivotal cellular activities. Histone crotonylation, as an essential reversible modification type, has become a research focus, laying the foundation for developing targeted therapeutic strategies (28).

This study utilized single-cell sequencing to resolve the impact of FOLFOX-Bev treatment on CC and further classified 8 cell types in light of Kcr modification levels. By systematically comparing communication characteristics between different groups, key findings were revealed. In the communication network of 9 cell types (including malignant cells classified into two types by Kcr level), T cells, fibroblasts, and endothelial cells demonstrated strong communication intensity and high frequency. Under the influence of high Kcr levels, L-R interactions involving T cells, fibroblasts, and endothelial cells were altered. In T cells, the CXCL16-CXCR6 and CXCL14-CXCR4 axes were inhibited, potentially impairing resident memory T cell (Trm) function and migration localization, thereby weakening immune responses (29,30). The HB-EGF/EGFR axis in the tumor microenvironment has been confirmed to be an important interaction center between tumor cells and fibroblasts (31). However, the EDN3-EDNRA axis was activated, and its function requires further study. In endothelial cells, the BMP2-BMPR2 axis was activated. A related study reported that it promotes angiogenesis in lung metastatic lesions in a mouse model of breast cancer (32). Based on these results, we speculate that higher Kcr levels enable specific ligands to bind receptors on T cells, fibroblasts, and endothelial cells, thereby regulating TF activities within these receptor cells and collectively shaping a unique microenvironmental state.

In the context of FOLFOX-Bev treatment, integrating TCGA data and utilizing the prognostic model, we identified 4 key genes: USP53, UACA, C7orf50, and IDH2. USP53 is a member of the ubiquitin-specific peptidase family, which is involved in numerous life processes (33). USP53 is a member of the ubiquitin-specific protease (USP) family, which is involved in numerous life activities. Given that the deubiquitination function of USPs can inherently affect protein stability, transcriptional regulation, cellular stress responses, and immune signaling pathways, USP53 may be involved in drug resistance and immune regulation (34). Furthermore, USP53 plays a role in DNA damage and the cell cycle, and is associated with drug resistance in lung cancer treatment and prognosis in CRC, with patients expressing low levels having a poorer prognosis (33). UACA encodes a protein containing ankyrin repeat sequences and coiled-coil domains, also known as Nucling, associated with apoptosis (35). In breast cancer patients taking Par-4-dependent chemotherapy drugs, high expression of the UACA gene was significantly associated with poorer overall survival, suggesting that UACA may be involved in the regulation of immune response and drug resistance mechanisms (36). However, the role of UACA in CRC remains unclear. Additionally, direct literature evidence elucidating its relationship with Kcr is currently scarce, necessitating further in-depth investigation. C7orf50 encodes cholesin hormone, which inhibits hepatic cholesterol synthesis, leading to reduced circulating cholesterol levels (37). Recent studies found that low cholesterol biosynthesis favors EMT maintenance and affects tumor molecular subtypes and disease-free survival in CC patients (38). C7orf50 may influence hepatic cholesterol synthesis by encoding the hormone cholesin, thereby participating in CC subtyping and patient survival regulation. However, investigations into the additional functions of C7orf50 remain sparse, particularly in the context of CRC. Therefore, our findings indicate that future inquiries should center on Kcr or drug resistance. IDH2, the mitochondrial isoform of isocitrate dehydrogenase, when mutated in multiple cancer types, leads to aberrant production of the oncometabolite D-2-hydroxyglutarate, profoundly affecting epigenetic, differentiation, and metabolic programs and making mutant IDH an ideal therapeutic target (39). IDH2 has been confirmed to be appreciably increased in CRC cells, actively promoting in vitro cell growth and in vivo tumor development (40). We also found IDH2 is one of the direct targets of Kcr modification regulated by Acox2, and its downregulated Kcr level may be a key molecular event in dysfunction, metabolic disorders, and hepatocarcinogenesis (41). Notably, although previous research has suggested that IDH2 may promote tumor growth in CRC, our multivariate Cox model showed a protective coefficient for IDH2 (HR <1). This difference is not contradictory. Prognostic models reflect conditional correlations in specific clinical cohorts after controlling for other variables, and are not equivalent to the direct oncogenic or tumor-suppressive causal effects of a single gene. IDH2 may exhibit different roles under different molecular subtypes, metabolic states, and treatment stresses. This study first discovered its high expression in FOLFOX-Bev drug-resistant subpopulations, potentially a pivotal target for Kcr-driven resistance. However, its specific mechanism and related pathways require experimental validation, suggesting future exploration into IDH2-mediated CC chemotherapy resistance and Kcr modification mechanisms.

Immune infiltration studies found increased pDCs infiltration in the high-risk group, accompanied by elevated HLA and MHC content. Concurrently, neutrophils and NK cells infiltration decreased, indicating the high-risk group has relatively stronger antigen-presenting capacity but lower involvement of innate immunity in disease risk control. Immune checkpoint analysis results revealed elevated PDCD1 expression and reduced CD40LG in the high-risk group. PDCD1, more commonly known as PD-1, is often considered to be part of the activation of anti-tumor immune responses (especially T cell responses) (42,43). CD40LG encodes CD40L on T cell surfaces (44). Research indicates CD40L is the core switch for DC transition from an immune-resting state to an immune-activated state, determining dendritic cell (DC) maturity and functional phenotype, indirectly shaping T cell response quality (effector function, memory formation), and ultimately affecting anti-tumor immune efficacy (45). As intervenable immune checkpoints, their association with protein expression and immunotherapy response requires further validation. Notably, the high-risk group in this study displayed a lower mutation burden and possessed a diminished capacity to forecast responses to immunotherapy. This phenomenon may reveal the complexity of the immune escape mechanism in CC. In a review study by Wang et al., Kcr can drive immunotherapy resistance by regulating gene expression to remodel the immune microenvironment, upregulating immune checkpoint molecules, and inhibiting co-stimulatory signals (46). Future research should validate the relationship between Kcr and immunotherapy response in larger cohorts and explore its potential as a biomarker or therapeutic target.

We then screened some potentially effective drugs for high-/low-risk group patients. In high-risk group predictions, Erlotinib_1168, AZD3759_1915, KRAS G12C Inhibitor 12_1855, and Sinularin_1838 have relevant research. Erlotinib_1168, an erlotinib (EGFR inhibitor) derivative, primarily demonstrates clinical efficacy in colitis-associated cancer models through potent tumor prevention and key signaling pathway regulation (47). AZD3759, an EGFR inhibitor, shows high efficacy and safety in the first-line treatment of non-small cell lung cancer (48). KRAS G12C Inhibitor is a compound under clinical investigation, suggesting pronounced preclinical anti-tumor activity in CC research (49). Sinularin has only been reported in prostate cancer cells, primarily inducing autophagy-dependent cell death by activating ULK1 and enhancing the FOXO3-ATG4A axis, but no detailed reports exist in CC (50). In the low-risk group, we selected Dasatinib_1079, Sepantronium.bromide_1941, AZD7762_1022, AZD1332_1463, IGF1R_3801_1738, and WEHI.1997 as candidate therapeutic drugs. Dasatinib, a BCR-ABL inhibitor and FDA-approved small molecule protease inhibitor, inhibits cell migration via focal adhesions and E-cadherin (51,52). Sepantronium bromide (also known as YM-155), a Survivin inhibitor, induces apoptosis and autophagy by downregulating Survivin, thereby inhibiting CC development (53). AZD7762 targets LXRα, LXRβ, and FXR pathways and is a promising anticancer candidate (54). However, studies regarding AZD133, IGF1R_3801_1738, WEHI.1997, TAF_5496_1732, and OF_1853 within CRC are currently absent. While direct evidence linking these agents to Kcr and drug resistance remains elusive, this study predicts that these candidate drugs may act on CRC-specific pathways based on the association analysis between Kcr-related characteristics and drug resistance phenotypes, potentially providing new strategies for overcoming drug resistance.

First, this study, based on transcriptomic and phenotypic analysis, primarily revealed the association between Kcr-related molecular states and FOLFOX-Bev treatment response, rather than directly proving a one-way causal effect of Kcr modification on drug resistance. Future research should systematically evaluate the direct relationship between Kcr levels and drug-resistant cells in larger single-cell cohorts to more rigorously examine the potential causal link. Second, based on a pre-defined set of Kcr genes, the AUCell score is an indirect inference indicator and cannot accurately reflect the actual protein modification level or site specificity. Direct validation using crotonylomics or immunohistochemistry techniques is required. Immune infiltration analysis relies on algorithmic prediction and lacks experimental evidence, such as flow cytometry or multicolor immunohistochemistry. Furthermore, the depth of the mechanism and the translational chain need to be strengthened. Our computational model-based screening of core genes and candidate drugs (such as erlotinib, AZD3759, dasatinib, and sepantronium bromide) requires functional validation in in vitro and in vivo experiments to enhance the reliability of our findings. In addition, this study, based on a multi-level bioinformatics analysis workflow, involves numerous parallel statistical tests, inevitably introducing the risk of false positives due to multiple hypothesis testing. Although conventional statistical threshold controls were employed in each analytical step, the possibility of cumulative statistical bias could not be completely eliminated. Simultaneously, the construction of the prognostic model relies on multi-step feature selection and regression analysis. While we mitigated the risk of overfitting to some extent through LASSO regression, independent validation cohorts, and multivariate Cox analysis, the model’s stability and generalization ability still require further validation in larger-scale, prospective cohorts.

Conclusions

In summary, this study is the first to reveal at single-cell resolution that Kcr modification drives FOLFOX-Bev resistance in CC by remodeling the immune microenvironment. The prognostic model constructed based on four key genes (USP53/UACA/C7orf50/IDH2) can accurately stratify patient risk and provide novel insights for targeted therapy and combined immunotherapy strategies in CC. Six candidate drug pairs were identified for high- and low-risk patients, offering potential strategies to overcome drug resistance.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-703/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-703/prf

Funding: This work was supported by Fujian Medical University QiHang Foundation General Project (grant No. 2024QH1330).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-703/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Morgan E, Arnold M, Gini A, et al. Global burden of colorectal cancer in 2020 and 2040: incidence and mortality estimates from GLOBOCAN. Gut 2023;72:338-44. [Crossref] [PubMed]
2. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
3. Schlechter BL. Management of Rectal Cancer. Hematol Oncol Clin North Am 2022;36:521-37. [Crossref] [PubMed]
4. Tonini V, Zanni M. Why is early detection of colon cancer still not possible in 2023? World J Gastroenterol 2024;30:211-24. [Crossref] [PubMed]
5. Yang M, Yang C, Ma D, et al. Single-cell analysis reveals cellular reprogramming in advanced colon cancer following FOLFOX-bevacizumab treatment. Front Oncol 2023;13:1219642. [Crossref] [PubMed]
6. Seufferlein T, Lausser L, Stein A, et al. Prediction of resistance to bevacizumab plus FOLFOX in metastatic colorectal cancer-Results of the prospective multicenter PERMAD trial. PLoS One 2024;19:e0304324. [Crossref] [PubMed]
7. Lim SH, Cho HJ, Kim KM, et al. Comprehensive molecular analysis to predict the efficacy of chemotherapy containing bevacizumab in patients with metastatic colorectal cancer. Oncol Res 2023;31:855-66. [Crossref] [PubMed]
8. Tan M, Luo H, Lee S, et al. Identification of 67 histone marks and histone lysine crotonylation as a new type of histone modification. Cell 2011;146:1016-28. [Crossref] [PubMed]
9. Han F, Shen W, Zhang X, et al. Protein crotonylation in cancer: mechanisms, functions, and therapeutic potential. Cell Biol Toxicol 2025;42:10. [Crossref] [PubMed]
10. Hou L, Chen YJ, Zhong Q, et al. Function and mechanism of lysine crotonylation in health and disease. QJM 2024;117:695-708. [Crossref] [PubMed]
11. Liao M, Sun X, Zheng W, et al. LINC00922 decoys SIRT3 to facilitate the metastasis of colorectal cancer through up-regulation the H3K27 crotonylation of ETS1 promoter. Mol Cancer 2023;22:163. [Crossref] [PubMed]
12. Mangiola S, Doyle MA, Papenfuss AT. Interfacing Seurat with the R tidy universe. Bioinformatics 2021;37:4100-7. [Crossref] [PubMed]
13. Aran D, Looney AP, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 2019;20:163-72. [Crossref] [PubMed]
14. Holland CH, Tanevski J, Perales-Patón J, et al. Robustness and applicability of transcription factor and pathway analysis tools on single-cell RNA-seq data. Genome Biol 2020;21:36. [Crossref] [PubMed]
15. Zhang Y, Liu T, Hu X, et al. CellCall: integrating paired ligand-receptor and transcription factor activities for cell-cell communication. Nucleic Acids Res 2021;49:8520-34. [Crossref] [PubMed]
16. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
17. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
18. Shi Y, Wang Y, Dong H, et al. Crosstalk of ferroptosis regulators and tumor immunity in pancreatic adenocarcinoma: novel perspective to mRNA vaccines and personalized immunotherapy. Apoptosis 2023;28:1423-35. [Crossref] [PubMed]
19. Xiao J, Wang J, Zhou C, et al. Development and Validation of a Propionate Metabolism-Related Gene Signature for Prognostic Prediction of Hepatocellular Carcinoma. J Hepatocell Carcinoma 2023;10:1673-87. [Crossref] [PubMed]
20. Zhao R, Wang W, Pan L, et al. The prognostic value and response to immunotherapy of immunogenic cell death-associated genes in breast cancer. Front Oncol 2023;13:1047973. [Crossref] [PubMed]
21. Liu S, Wang Z, Zhu R, et al. Three Differential Expression Analysis Methods for RNA Sequencing: limma, EdgeR, DESeq2. J Vis Exp 2021; [PubMed]
22. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
23. Wang L, Wang D, Yang L, et al. Cuproptosis related genes associated with Jab1 shapes tumor microenvironment and pharmacological profile in nasopharyngeal carcinoma. Front Immunol 2022;13:989286. [Crossref] [PubMed]
24. Zhu S, Cheng Q, Zou M, et al. Combining bulk and scRNA-seq to explore the molecular mechanisms governing the distinct efferocytosis activities of a macrophage subpopulation in PDAC. J Cell Mol Med 2024;28:e18266. [Crossref] [PubMed]
25. 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:bbab260. [Crossref] [PubMed]
26. Nwabo Kamdje AH, Dongmo Fogang HP, Mimche PN. Role of epigenetic in cancer biology, in hematologic malignancies and in anticancer therapy. Front Mol Med 2024;4:1426454. [Crossref] [PubMed]
27. Mahmoodi Chalbatani G, Gharagouzloo E, Malekraeisi MA, et al. The integrative multi-omics approach identifies the novel competing endogenous RNA (ceRNA) network in colorectal cancer. Sci Rep 2023;13:19454. [Crossref] [PubMed]
28. Xie JY, Ju J, Zhou P, et al. The mechanisms, regulations, and functions of histone lysine crotonylation. Cell Death Discov 2024;10:66. [Crossref] [PubMed]
29. Bao N, Fu B, Zhong X, et al. Role of the CXCR6/CXCL16 axis in autoimmune diseases. Int Immunopharmacol 2023;121:110530. [Crossref] [PubMed]
30. Salve BG, Sharma S, Vijay N. Evolutionary diversity of CXCL16-CXCR6: Convergent substitutions and recurrent gene loss in sauropsids. Immunogenetics 2024;76:397-415. [Crossref] [PubMed]
31. Giangreco G, Rullan A, Naito Y, et al. Cancer cell - Fibroblast crosstalk via HB-EGF, EGFR, and MAPK signaling promotes the expression of macrophage chemo-attractants in squamous cell carcinoma. iScience 2024;27:110635. [Crossref] [PubMed]
32. Harry JA, Skebo SI, Cole DV, et al. Bmpr2 Loss in the Vascular Endothelium Enhances Neoangiogenesis and Growth of Lung Metastatic Lesions in a Mouse Model of Breast Cancer. J Cell Biochem 2025;126:e70071. [Crossref] [PubMed]
33. Xia G, Guo Y, Zhang J, et al. An Overview of the Deubiquitinase USP53: A Promising Diagnostic Marker and Therapeutic Target. Curr Protein Pept Sci 2024;25:708-18. [Crossref] [PubMed]
34. Gao H, Xi Z, Dai J, et al. Drug resistance mechanisms and treatment strategies mediated by Ubiquitin-Specific Proteases (USPs) in cancers: new directions and therapeutic options. Mol Cancer 2024;23:88. [Crossref] [PubMed]
35. Matsui T, Sakamaki Y, Nakashima S, et al. Rab39 and its effector UACA regulate basolateral exosome release from polarized epithelial cells. Cell Rep 2022;39:110875. [Crossref] [PubMed]
36. Zhu Q, Schultz E, Long J, et al. UACA locus is associated with breast cancer chemoresistance and survival. NPJ Breast Cancer 2022;8:39. [Crossref] [PubMed]
37. Hu X, Chen F, Jia L, et al. A gut-derived hormone regulates cholesterol metabolism. Cell 2024;187:1685-1700.e18. [Crossref] [PubMed]
38. Aulas A, Liberatoscioli ML, Finetti P, et al. Low cholesterol biosynthesis favors epithelial-to-mesenchymal transition maintenance and influences tumor molecular subtyping and disease-free survival in colon cancer patients. Cancer Commun (Lond) 2022;42:793-7. [Crossref] [PubMed]
39. Pirozzi CJ, Yan H. The implications of IDH mutations for cancer development and therapy. Nat Rev Clin Oncol 2021;18:645-61. [Crossref] [PubMed]
40. White K, Someya S. The roles of NADPH and isocitrate dehydrogenase in cochlear mitochondrial antioxidant defense and aging. Hear Res 2023;427:108659. [Crossref] [PubMed]
41. Zhang Y, Chen Y, Zhang Z, et al. Acox2 is a regulator of lysine crotonylation that mediates hepatic metabolic homeostasis in mice. Cell Death Dis 2022;13:279. [Crossref] [PubMed]
42. Lin X, Kang K, Chen P, et al. Regulatory mechanisms of PD-1/PD-L1 in cancers. Mol Cancer 2024;23:108. [Crossref] [PubMed]
43. Ibrahiem AT, Eladl E, Toraih EA, et al. Prognostic Value of BRAF, Programmed Cell Death 1 (PD1), and PD Ligand 1 (PDL1) Protein Expression in Colon Adenocarcinoma. Diagnostics (Basel) 2023;13:237. [Crossref] [PubMed]
44. Santaniemi W, Pernaa N, Glumoff V, et al. CD40LG Triplication Associates with Immune Dysregulation and Exhaustion. J Clin Immunol 2023;43:323-6. [Crossref] [PubMed]
45. de Mey W, Locy H, De Ridder K, et al. An mRNA mix redirects dendritic cells towards an antiviral program, inducing anticancer cytotoxic stem cell and central memory CD8(+) T cells. Front Immunol 2023;14:1111523. [Crossref] [PubMed]
46. Wang X, Qu Y, Li Z, et al. Histone crotonylation in tumors Mol Clin Oncol 2025;22:39. (Review). [Crossref] [PubMed]
47. Liu M, Zhong XS, Krishnachaitanya SS, et al. Erlotinib suppresses tumorigenesis in a mouse model of colitis-associated cancer. Biomed Pharmacother 2024;175:116580. [Crossref] [PubMed]
48. Maggie Liu SY, Dong XR, Wang Z, et al. Efficacy, safety and dose selection of AZD3759 in patients with untreated EGFR-mutated non-small-cell lung cancer and central nervous system metastases in China (CTONG1702-Arm 8): a multi-center, single-arm, phase 2 trial. EClinicalMedicine 2023;64:102238. [Crossref] [PubMed]
49. Tang D, Kang R. Glimmers of hope for targeting oncogenic KRAS-G12D. Cancer Gene Ther 2023;30:391-3. [PubMed]
50. Meng XY, Li Y, Yan ZJ, et al. Sinularin induces autophagy-dependent cell death by activating ULK1 and enhancing FOXO3-ATG4A axis in prostate cancer cells. Sci Rep 2025;15:15875. [Crossref] [PubMed]
51. Lu YW, Hou XL, Koo HM, et al. Dasatinib suppresses collective cell migration through the coordination of focal adhesion and E-cadherin in colon cancer cells. Heliyon 2024;10:e23501. [Crossref] [PubMed]
52. Roskoski R Jr. Properties of FDA-approved small molecule protein kinase inhibitors: A 2024 update. Pharmacol Res 2024;200:107059. [Crossref] [PubMed]
53. Liu N, Li Y, Luo G, et al. SIRT6 suppresses colon cancer growth by inducing apoptosis and autophagy through transcriptionally down-regulating Survivin. Mitochondrion 2024;78:101932. [Crossref] [PubMed]
54. Singha B, Gogoi PP, Longkumer P, et al. A holistic computational exploration of AZD7762 as a potent selective modulator of LXRα, LXRβ and FXR: An underexplored pathway in cancer therapeutics. Comput Biol Med 2025;194:110433. [Crossref] [PubMed]