A radiotherapy resistance-related prognostic signature predicts survival and the immune landscape in rectal cancer

Introduction

Colorectal cancer (CRC) is the third most common cancer and the second leading cause of cancer-related death worldwide, and thus represents a major health issue (1). Rectal cancer (RC) accounts for nearly 30% of CRC cases. Approximately one-third of RC patients are diagnosed with locally advanced rectal cancer (LARC) (2).

Neoadjuvant radiotherapy, combined with surgery and chemotherapy, has significantly improved local control and sphincter preservation in patients with RC. The introduction of total neoadjuvant therapy (TNT) (3), which delivers all systemic chemotherapy before surgery, has further optimized treatment regimens for LARC. Multiple randomized trials, such as the RAPIDO and PRODIGE 23 studies, have demonstrated that TNT significantly improves pathological complete response rates (up to 28–30%) and disease-free survival compared to standard neoadjuvant chemoradiotherapy (3,4).

Despite ongoing innovations in multimodal treatment regimens over the past few decades, nearly 30% of patients with LARC develop intrinsic or acquired resistance to radiotherapy (6). This resistance leads to tumor locoregional recurrence and distant relapse, ultimately contributing to a poor prognosis. Identifying reliable biomarkers for predicting radiotherapy response and long-term survival remains a major clinical challenge.

Evidence suggests that radiotherapy resistance in RC is a multifactorial process involving both tumor-intrinsic and tumor-extrinsic mechanisms (6,7). Enhanced DNA damage repair capacity allows tumor cells to tolerate radiation-induced genomic stress, thereby attenuating treatment efficacy (9). In parallel, epithelial-mesenchymal transition (EMT) and increased cellular plasticity have been shown to promote radioresistant phenotypes and facilitate tumor invasion and recurrence (10).

Beyond tumor-intrinsic alterations, the tumor microenvironment plays a critical role in shaping radiotherapy response. Aberrant cytokine signaling, chronic inflammation, and stromal activation can suppress radiation-induced antitumor immunity, leading to immune exclusion and therapeutic failure (11). Increasing evidence also implicates angiogenesis, hypoxia, and immune checkpoint (CP)-mediated immunosuppression in adaptive resistance following radiotherapy (12-14). These complex interactions highlight the significance of immune remodeling and cellular heterogeneity in determining radiotherapy response, underscoring the need for further investigation into the mechanisms underlying radiotherapy resistance.

Despite these advances, clinically applicable biomarkers that integrate radiotherapy resistance mechanisms with prognostic prediction have yet to be established. Moreover, the cellular heterogeneity and immune microenvironment underlying radiotherapy resistance at single-cell resolution remain incompletely understood. Accordingly, comprehensive integrative analyses are needed to clarify the links between radiotherapy resistance, prognosis, and tumor microenvironment characteristics, and to identify effective predictive markers.

In this study, we performed integrative transcriptomic and single-cell bioinformatics analyses to identify radiotherapy resistance-radiotherapy resistance-related genes and to develop a radiotherapy-related prognostic signature for RC. We further examined the functional pathways, immune landscape, and cellular heterogeneity associated with the signature genes, providing insights into their association with radiotherapy resistance and prognostic assessment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0554/rc).

Methods

Data acquisition and preprocessing

RC transcriptomic data [The Cancer Genome Atlas Rectum Adenocarcinoma (TCGA-READ)] were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) in transcripts per million (TPM) format. The study population comprised patients with RC with available transcriptomic and clinical follow-up data. Patients were included in the study if they had complete gene expression profiles and survival information. Patients were excluded from the analysis if key clinical variables or follow-up data were unavailable. A total of 154 patients with complete clinical and follow-up information were included as the training cohort. Gene expression data were normalized prior to analysis, and no additional imputation for missing values was performed.

The radiotherapy-related RC dataset GSE35452 was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE35452) using the GEO query R package (15). The GSE35452 (GPL570 platform) cohort included 46 patients, comprising 24 responders and 22 non-responders to preoperative radiotherapy. The GSE87211 dataset (GPL13497) was also obtained from the GEO and used as an independent validation cohort. The GEO datasets were normalized and standardized using the limma R package. Single-cell RNA sequencing (scRNA-seq) data were obtained from the GEO dataset GSE144735, which included six CRC tumor samples, for subsequent single-cell analysis.

The datasets used in this study are publicly available, and the accrual and follow-up periods were defined according to the original studies from which the data were obtained. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Differential expression analysis

Differentially expressed genes (DEGs) associated with radiotherapy response were identified by comparing the radiotherapy-responsive and non-responsive groups in the GSE35452 dataset using the limma R package (16). Genes with a |log₂ fold change| ≥0.585 and a P value <0.05 were considered statistically significant.

Weighted gene correlation network analysis (WGCNA)

A WGCNA was performed to identify the gene modules associated with the radiotherapy-related phenotypes (17). The soft-thresholding power was determined using the pickSoftThreshold function to construct a scale-free network. The adjacency matrix was transformed into a topological overlap matrix, followed by hierarchical clustering to identify co-expression modules. Modules with an absolute correlation coefficient >0.4 and a P value <0.05 were considered significantly associated with clinical phenotypes.

Construction of the prognostic model

Candidate predictors included radiotherapy response-related genes identified through differential expression analysis and WGCNA, and their expression levels were treated as continuous variables in the subsequent analyses. A univariate Cox proportional hazards regression analysis was first performed to examine the association between radiotherapy-related genes and overall survival (OS) in TCGA-READ cohort, with P<0.05 as the selection criterion. Least absolute shrinkage and selection operator (LASSO) regression was then applied using the glmnet package (18) to eliminate multicollinearity and select key variables (set.seed =1, family = “cox”). Subsequently, a multivariate Cox regression analysis with stepwise selection was conducted to identify the independent prognostic genes and construct the final prognostic model. The risk score was calculated as follows:

Riskscore=∑i=1ncoefi∗expi[1]

Patients were stratified into high- and low-risk groups based on the calculated risk scores. Kaplan-Meier survival analysis and log-rank tests were used to compare OS between the groups. Time-dependent receiver operating characteristic (ROC) curves were generated using the timeROC package (19) to evaluate the predictive accuracy of the model. The GSE87211 dataset was used as an external validation cohort.

Functional enrichment analysis

A Gene Ontology (GO) enrichment analysis, including biological process (BP), molecular function (MF), and cellular component (CC), and a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed on the genes associated with the radiotherapy resistance risk score using the clusterProfiler R package. An adjusted P value <0.05 was considered statistically significant. A gene set enrichment analysis (GSEA) (20) was conducted using the Hallmark gene set (“h.all.v7.0.entrez.gmt”) and the KEGG gene set (“c2.cp.kegg.v7.0.entrez.gmt”) downloaded from the Molecular Signatures Database (MSigDB), with a P value <0.05 indicating significant enrichment.

Immune-related analysis

Immune cell infiltration in TCGA-READ cohort was assessed using the IOBR R package (21), which integrates the CIBERSORT (22), EPIC (23), and ESTIMATE (24) algorithms. Differences in immune cell composition between the high- and low-risk groups were analyzed. Further, the immunophenoscore (IPS) was computed to evaluate tumor immunogenicity based on the major histocompatibility complex (MHC), immune CPs, effector cells (ECs), suppressor cells (SCs), and the composite score (average Z score).

Somatic mutation and immune CP analysis

Somatic mutation data from TCGA-READ cohort were downloaded and processed to evaluate the mutational landscape between the high- and low-risk groups. Mutation annotation format (MAF) files were analyzed to summarize overall mutation frequency, variant classification, and mutation spectrum using standard bioinformatics approaches. A differential mutation analysis between subgroups was performed to identify group-specific mutated genes. Immune CP gene expression levels, including CD274, PDCD1LG2, PDCD1, LAG3, CD70, and CTLA4, were extracted from the transcriptomic data and compared between the two risk groups.

Nomogram construction

A prognostic nomogram was constructed based on a multivariate Cox regression analysis incorporating risk score, node (N) stage, and age in TCGA-READ cohort. The nomogram was used to predict patient prognosis, and its performance was evaluated using calibration curves.

scRNA-seq analysis

The Seurat R package (25) was used to process scRNA-seq data obtained from GSE144735. Cells with mitochondrial gene content >30%, fewer than 500 detected genes, or more than 5,000 detected genes were eliminated. After normalization, 2,000 highly variable genes were identified and a principal component analysis was performed. Thirty principal components were selected for uniform manifold approximation and projection (UMAP) dimensionality reduction. Cells were subsequently clustered using the FindClusters function at a resolution of 0.8, resulting in the identification of 19 distinct clusters. Cell types were annotated based on canonical marker genes.

Cell-cell communication analysis

The cell-cell communication analysis was performed using the CellChat R package (26) (http://www.cellchat.org/) to infer and quantify intercellular interactions based on known ligand-receptor pairs. Interaction probability and permutation tests were applied to identify significant interactions. The number and intensity of significant ligand-receptor interactions between cell types were visualized using heatmaps and interaction graphs.

Pseudotime trajectory analysis

Tumor epithelial cells were identified from the dataset GSE144735. A pseudotime analysis was then performed using the R package Monocle (27) to infer a branched differentiation trajectory. Gene expression dynamics were assessed across the pseudotime axis.

Statistical analysis

The primary outcome of this study was OS, defined as the time from diagnosis to death from any cause or last follow-up. Gene expression values were analyzed as continuous variables without categorization unless otherwise specified (e.g., risk group stratification based on the median risk score). Continuous variables between two groups were compared using the Mann-Whitney U test. Survival analyses were performed using the survival and survminer R packages. Univariate and multivariate Cox regression analyses were performed to identify prognostic factors. All statistical analyses were conducted using R software (version 4.2.3). All statistical tests were two-sided, and a P value <0.05 was considered statistically significant.

Results

Identification of radiotherapy response-related gene modules in RC

A WGCNA was performed on the GSE35452 dataset to identify the gene modules associated with radiotherapy response in RC. An appropriate soft-thresholding power was chosen to approximate scale-free topology (Figure 1A,1B), and the genes were subsequently clustered into distinct co-expression modules using hierarchical clustering and dynamic tree cutting (Figure 1C).

Figure 1 WGCNA and differential expression analysis of radiotherapy response-related genes. (A-C) WGCNA of the GSE35452 dataset was performed, and genes were clustered into distinct modules as shown by the gene dendrogram. (D) Module-trait correlation heatmap for radiotherapy response. (E) Relationship between module membership and gene significance in the cyan module. (F) Volcano plot of DEGs between radiotherapy responders and non-responders. (G) Heatmap showing expression patterns of the identified DEGs across samples. DEGs, differentially expressed genes; FC, fold change; WGCNA, weighted gene correlation network analysis.

A module-trait relationship analysis indicated that the cyan module showed the strongest association with radiotherapy response status (cor =0.42, P=0.005; Figure 1D). The cyan module contained 487 genes, and within this module, module membership was strongly correlated with gene significance for radiotherapy response (cor =0.53, P=1.3×10⁻³⁶), suggesting that the cyan module genes were more likely to be biologically relevant to treatment response (Figure 1E).

A differential expression analysis between the radiation responders and non-responders identified 61 DEGs meeting the criteria of |log₂ fold change| ≥0.585 and P<0.05. Within this dataset, the radiation-resistant group exhibited 57 up-regulated genes and four down-regulated genes (Figure 1F). A heatmap of these 61 DEGs revealed significant differential expression patterns between the two groups (Figure 1G), confirming the presence of radiation response-related transcriptional characteristics. By integrating the radiation-related DEGs with the genes from the WGCNA modules, we ultimately obtained and screened candidate genes based on survival characteristics.

Construction and validation of a radiotherapy-related prognostic signature

In TCGA-READ cohort, univariate Cox regression identified 25 radiotherapy-related genes significantly associated with OS (P<0.05; Figure 2A). Most candidates exerted protective effects [hazard ratio (HR) <1]; however, PSCA presented as a risk factor [HR =1.354, 95% confidence interval (CI): 1.079–1.699, P=0.009].

Figure 2 Construction of a radiotherapy-related prognostic signature. (A) Forest plot of univariate Cox regression identifying radiotherapy-related genes associated with overall survival in TCGA-READ cohort. (B,C) LASSO Cox regression analysis for feature selection. (D,E) Multivariate Cox regression analysis and corresponding coefficients of the four-gene prognostic model. (F-H) Risk score distribution, patient stratification, survival status, and expression patterns of the four genes between risk groups. (I,J) Kaplan-Meier survival curve and time-dependent ROC analysis in TCGA-READ cohort. (K,L) Kaplan-Meier survival curve and ROC analysis in the external validation cohort GSE87211. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA-READ, The Cancer Genome Atlas Rectum Adenocarcinoma.

These prognostic candidates were further subjected to LASSO Cox regression (Figure 2B,2C), followed by multivariate Cox regression to construct the final model. A four-gene signature comprising CUTA, IZUMO2, PALB2, and PSCA was established (Figure 2D). In the multivariate model, IZUMO2 (HR =0.572, P=0.03) and PALB2 (HR =0.553, P=0.001) remained independent protective factors, while PSCA (HR =1.400, P=0.03) remained an independent risk factor; CUTA exerted a protective effect but this effect was not statistically significant (HR =0.769, P=0.16) (Figure 2D). Consistently, the model assigned a positive coefficient to PSCA and negative coefficients to CUTA, IZUMO2 and PALB2 (Figure 2E).

Based on the multivariate Cox coefficients, a risk score was calculated for each patient as a linear combination of the expression levels of the four genes. Patients were stratified into high- and low-risk groups according to the median value of the calculated risk score (Figure 2F). Patients with higher risk scores exhibited a higher proportion of deaths and shorter survival times (Figure 2G), and the four model genes displayed distinct expression patterns between the risk groups (Figure 2H).

The Kaplan-Meier analysis revealed that OS was significantly worse in the high-risk group in TCGA-READ training cohort (P<0.001; Figure 2I). Time-dependent ROC analysis yielded area under the curve (AUC) values of 0.724, 0.760, and 0.842 at 1, 2, and 3 years, respectively (Figure 2J). Notably, the signature showed consistent prognostic performance in the external validation cohort GSE87211, where the high-risk group exhibited worse OS (P=0.002; Figure 2K), with AUC values of 0.765, 0.797, and 0.744 at 1, 2, and 3 years, respectively (Figure 2L).

Functional characteristics associated with the radiotherapy-related risk score

To elucidate the biological mechanisms underlying the radiotherapy resistance-related risk score, functional enrichment analyses were performed on the risk score-related genes. The GO analysis demonstrated significant enrichment across the BP, CC, and MF categories, indicating that the risk score reflects coordinated alterations in multiple functional categories rather than isolated molecular events (Figure 3A).

Figure 3 Functional enrichment analysis of risk score-associated genes. (A) GO enrichment analysis across the BP, CC, and MF categories. (B) KEGG pathway enrichment analysis of the risk score-associated genes. (C-G) KEGG-based GSEA comparing high- and low-risk groups. (H-L) Hallmark GSEA showing pathways differentially enriched between the high- and low-risk groups. BP, biological process; CC, cellular component; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

The KEGG pathway analysis further showed that the risk score-related genes were enriched in pathways related to signal transduction and cellular dynamics, including mitogen-activated protein kinase (MAPK) signaling, Ras/Ras-related protein 1 (Rap1) signaling, calcium signaling, endocytosis, focal adhesion, regulation of the actin cytoskeleton, as well as pathways involved in cell cycle regulation and chemokine signaling (Figure 3B).

The GSEA using the KEGG gene sets revealed that the inflammation- and immune-related pathways were enriched in the high-risk group (Figure 3C,3F,3G), whereas the mismatch repair and DNA replication related pathways were enriched in the low-risk group (Figure 3C-3E).

Consistently, the Hallmark GSEA revealed enrichment of the EMT, inflammatory response, and angiogenesis gene sets in the high-risk group (Figure 3H,3J-3L). In contrast, MYC TARGETS V2 gene sets were enriched in the low-risk group (Figure 3H,3I). Overall, differences in immune/inflammatory and genome maintenance-related transcriptional patterns were observed between the high- and low-risk groups.

Immune microenvironment features associated with the risk score

Immune cell infiltration patterns across risk subgroups

To further characterize the immune landscape associated with the radiotherapy resistance-related risk score, immune and stromal features were compared between the high- and low-risk groups.

The ESTIMATE analysis showed that the high-risk group had higher ESTIMATE scores (P=0.01) and immune scores (P=0.004) than the low-risk group, whereas no significant difference was observed in the stromal scores between the groups (P=0.09) (Figure 4A-4C). Similarly, the integrative heatmap analysis based on CIBERSORT, EPIC, and ESTIMATE revealed distinct immune infiltration patterns and microenvironmental features between the two risk groups (Figure 4D). The IPS analysis further indicated differences in immunogenicity-related components. Compared with the high-risk group, the low-risk group exhibited significantly higher SC-IPS, CP-IPS, and AZ-IPS (Figure 4E).

Figure 4 Immune-related analyses stratified by the risk score. (A-C) Comparison of the ESTIMATE, immune, and stromal scores between the high- and low-risk groups. (D) Heatmap summarizing immune infiltration patterns estimated by CIBERSORT, EPIC, and ESTIMATE. (E) Comparison of IPS components, including the MHC, ECs, SCs, CPs, and AZ scores, between the high- and low-risk groups. (F) Differential expression of selected immune-related and epigenetic regulatory genes between the high- and low-risk groups. (G) Correlation analysis between risk score and immune regulatory genes. AZ; CPs, checkpoints; ECs, effector cells; IPS, immunophenoscore; MHC, major histocompatibility complex; SC, suppressor cells.

The expression analysis of the immune-related genes showed differential expression of several immunoregulatory and epigenetic regulators, including PDCD1, EZH2 and DNMT1 (Figure 4F). In addition, the correlation analysis revealed associations between the risk score and multiple immune CP- and immune regulation-related genes (Spearman correlation; Figure 4G).

Overall, these results indicated that high-risk tumors show higher inferred immune infiltration, while low-risk tumors display higher immunogenicity-related scores and distinct immune regulatory profiles, which may be relevant to treatment response heterogeneity.

Genomic mutation landscape and immune CP profiling of risk subgroups

To further elucidate the genomic and immune characteristics underlying the risk stratification, somatic mutation and immune CP analyses were performed in TCGA-READ cohort. Somatic mutation profiling was conducted in the high-risk (n=62, Figure 5A) and low-risk (n=63, Figure 5B) subgroups. Genomic alterations were detected in 85.48% and 93.65% of samples, respectively. The most frequently mutated genes were APC and TP53, followed by TTN and KRAS, with missense mutation being the predominant alteration type, accompanied by nonsense mutations, frameshift insertions/deletions, and splice site variants.

Figure 5 Mutational landscape and immune checkpoint profiles between the high- and low-risk groups. (A,B) Top 10 most frequently mutated genes in the high- (A) and low-risk (B) groups. (C) Forest plot showing genes with significantly different mutation frequencies between the high- and low-risk groups. (D) Expression levels of representative immune checkpoint genes (CD274, PDCD1LG2, PDCD1, LAG3, CD70, and CTLA4) between the high- and low-risk groups. *; **, OR, odds ratio.

The differential mutation analysis between the risk groups revealed distinct mutational patterns (Figure 5C). The low-risk group was enriched in PIK3CA and SYNE1 mutations, while the high-risk group was characterized by mutations in CACNG3, SOGA1, ZNF521, PCLO, and ROS1, with several genes showing group-specific occurrence [odds ratio = infinity (Inf)].

In addition, the immune CP analysis demonstrated generally higher expression levels of immune-related genes in the high-risk group. PDCD1 (PD-1) was significantly upregulated, and the expression of LAG3 and CTLA4 was also increased (Figure 5D). We also systematically analyzed differences in tumor mutational burden (TMB) and microsatellite instability (MSI) status; however, no statistically significant differences were observed between the two cohorts (data not provided).

Independent prognostic value of the radiotherapy-related risk score

To assess whether the radiotherapy resistance-related risk score represents an independent prognostic factor, univariate and multivariate Cox regression analyses were performed incorporating clinicopathological variables. In the univariate analysis, risk score, N stage, metastasis (M) stage, and age were significantly associated with OS, whereas tumor (T) stage and gender were not (Figure 6A). Notably, the risk score was significantly associated with OS (HR =2.801; 95% CI: 1.827–4.294; P<0.001). Multivariate Cox regression confirmed that the risk score remained an independent prognostic factor (HR =2.532; 95% CI: 1.565–4.099; P<0.001), together with N stage (HR =3.051, P=0.03) and age (HR =7.974, P=0.045), while M stage lost statistical significance in the multivariate model (Figure 6B).

Figure 6 Nomogram incorporating the radiotherapy-related risk score. (A) Univariate Cox regression analysis of the clinicopathological variables and risk score. (B) Multivariate Cox regression identifying independent prognostic factors. (C) A nomogram was constructed based on the risk-score and clinical parameters. (D-F) Calibration curves of the nomogram for predicting OS at 1-, 2-, and 3-year in TCGA READ dataset. CI, confidence interval; HR, hazard ratio; OS, overall survival; TCGA-READ, The Cancer Genome Atlas Rectum Adenocarcinoma.

Based on these independent predictors, a prognostic nomogram integrating risk score, age, and N stage was constructed to estimate 1-, 2-, and 3-year OS probabilities (Figure 6C). The calibration curve analysis demonstrated good agreement between nomogram-predicted and observed OS at 1, 2, and 3 years (Figure 6D-6F).

Collectively, these results revealed that the radiotherapy resistance-related risk score remained independently associated with OS after adjustment for clinicopathological variables.

Survival analysis of prognostic signature genes and pathway enrichment associated with PSCA expression

The chromosomal locations of the four genes (CUTA, IZUMO2, PALB2, and PSCA) included in the prognostic signature are shown in Figure 7A. Kaplan-Meier survival analyses demonstrated that the expression levels of the individual signature genes were significantly associated with OS. Higher expression of CUTA was associated with improved OS (P=0.002; Figure 7B). Similarly, elevated expression of IZUMO2 and PALB2 was correlated with favorable OS (both P<0.001; Figure 7C,7D). In contrast, higher expression of PSCA was associated with poorer OS (P=0.03; Figure 7E).

Figure 7 Survival analyses and pathway enrichment of the four-gene signature. (A) Chromosomal locations of CUTA, IZUMO2, PALB2, and PSCA. (B-E) Kaplan-Meier OS analyses stratified by expression levels of individual signature genes. (F,G) DSS and PFI analyses based on PSCA expression. (H) Summary plot of Hallmark GSEA for selected pathways. (I,J) Representative enrichment plots for selected Hallmark pathways. DSS, disease-specific survival; GSEA, gene set enrichment analysis; OS, overall survival; PFI, progression-free interval.

Other survival endpoint data further confirmed the correlation between PSCA expression and poor prognosis. Patients with higher PSCA expression experienced worse disease-specific survival (P=0.02; Figure 7F) and a shorter progression-free interval (P=0.002; Figure 7G), indicating a consistent association between PSCA expression and unfavorable outcomes across these survival indicators.

A GSEA was performed to investigate the Hallmark pathways associated with PSCA expression levels (Figure 7H). The HALLMARK_KRAS_SIGNALING_DN gene set was enriched in the group with high PSCA expression (Figure 7I), whereas the HALLMARK_MYC_TARGETS_V2 gene set was enriched in the group with low PSCA expression (Figure 7J).

Single-cell analysis revealed PSCA+ epithelial cells in the tumor microenvironment

The scRNA-seq data from the GSE144735 dataset underwent UMAP dimensionality reduction analysis, resulting in the identification of 19 cell clusters with distinct transcriptional profiles (Figure 8A). Cells from various tumor samples exhibited a broad distribution and relative uniformity among the clusters, without obvious sample-specific segregation (Figure 8B).

Figure 8 Single-cell landscape and cell type-specific expression of prognostic genes. (A) UMAP plot of cell clusters from the GSE144735 dataset. (B) UMAP plot of cell distribution from different samples. (C) Cell type annotation based on canonical marker genes. (D) Dot plot showing the marker gene expression for each cell subpopulation. (E) Pie chart presenting the relative proportions of annotated cell types. (F) Heatmap of marker gene expression with associated GO-BP annotations for each cell subpopulation. (G-J) Feature plots showing cell type-specific expression of CUTA, IZUMO2, PALB2, and PSCA. BP, biological process; GO, Gene Ontology; UMAP, uniform manifold approximation and projection.

Based on the expression of marker genes, these cell clusters were classified into 10 major cell types, including epithelial cells, T or natural killer cells, fibroblasts, myeloid cells, endothelial cells, B cells, smooth muscle cells, plasma cells, mast cells, and glial cells (Figure 8C,8D). The relative proportions of each cell population are shown in Figure 8E. An integrated heatmap visualizing the results of the GO-BP analysis based on 10 different cell types revealed complex and diverse functions of the corresponding gene sets (Figure 8F).

UMAP dimensionality reduction analysis revealed cell type-specific expression patterns of the four prognostic genes (Figure 8G-8J). PSCA expression was predominantly observed in epithelial cells, whereas CUTA, IZUMO2 and PALB2 displayed more diffuse or relatively low expression across multiple cell types. These findings indicated that the prognostic information reflected by the four-gene signature may originate from distinct cell types in the tumor microenvironment, rather than from a single cell type alone.

Cell-cell communication and developmental trajectories of PSCA+ epithelial cells

To further investigate the function of the PSCA gene in tumor epithelial cells, based on the single-cell dataset GSE144735, the cells were classified as PSCA positive (PSCA+) or PSCA negative (PSCA−) according to the PSCA expression levels. A cell communication analysis was conducted, showing the quantity and intensity of cell communication between these groups and other cell types (Figure 9A,9B). The intensity of signaling sending and receiving in the PSCA+ epithelial cells exceeded that observed in the PSCA− epithelial cells (Figure 9C). The interactions of PSCA+ or PSCA− epithelial cells as signal senders with surrounding stromal cells were inferred.

Figure 9 Inferred epithelial cell communication analyses and pseudotime trajectories. (A,B) Overview of inferred cell-cell communication networks involving PSCA+ and PSCA− epithelial cells. (C) Comparison of outgoing and incoming interaction strengths between PSCA+ and PSCA− epithelial cells. (D-F) Ligand-receptor interactions of PSCA+ and PSCA− cells acting as signal emitters with immune and stromal cells. (G-I) Ligand-receptor interactions of PSCA+ and PSCA− cells acting as signal receivers with immune and stromal cells. (J,K) Distribution of PSCA expression at single-cell resolution. (L-N) Pseudotime trajectories illustrating inferred transcriptional heterogeneity between PSCA+ and PSCA− epithelial cells. UMAP, uniform manifold approximation and projection.

The bubble plot revealed that signaling pathways such as secreted phosphoprotein 1 (SPP1)-[integrin subunit alpha V (ITGAV) + integrin subunit beta 1 (ITGB1)] and macrophage migration inhibitory factor (MIF)-[CD74 + C-X-C motif chemokine receptor 4 (CXCR4)] (Figure 9D), midkine (MK)-(ITGA6 + ITGB1), and MK-nucleolin (NCL) (Figure 9E) were significantly enriched in the PSCA+ epithelial cells compared to the PSCA− epithelial cells. The heat map further revealed that the PSCA+ epithelial cells send higher signal intensities of MK and MIF-related molecules compared to the PSCA− cells in outgoing signaling patterns (Figure 9F).

The role of PSCA+ and PSCA− epithelial cells as signal receivers in cellular communication was investigated. Signaling pathways such as SPP1-(ITGAV + ITGB1), SPP1-CD44 (Figure 9G), MK-syndecan1 (SDC1), and MK-(ITGA6 + ITGB1) (Figure 9H) were more enriched in the PSCA+ epithelial cells compared to the PSCA− epithelial cells. The heatmap revealed that the PSCA+ epithelial cells received stronger signals of MK and insulin-like growth factor (IGF) compared to the PSCA− epithelial cells (Figure 9I).

At single-cell resolution, PSCA expression was confined to a subset of epithelial cells occupying distinct regions of transcriptional space (Figure 9J,9K). The pseudotime analysis demonstrated that the PSCA+ epithelial cells were preferentially enriched at higher pseudotime values, whereas the PSCA− epithelial cells were distributed across earlier and intermediate states (Figure 9L-9N).

Discussion

The radiotherapy response and prognosis of patients with RC are influenced by complex molecular and microenvironmental factors rather than isolated gene-level alterations. In this study, we integrated network-based transcriptomic analysis, survival modeling, immune profiling, and single-cell approaches to develop a radiotherapy resistance-related transcriptional signature and to explore its prognostic relevance.

By conducting a WGCNA on a radiotherapy-treated cohort, we identified a co-expression module connected with treatment response. After integrating the module genes with the DEGs, we obtained a candidate gene pool for the downstream prognostic model, thereby reducing the dependence on the effect of a single gene. From these candidates, we ultimately developed a four-gene prognostic signature comprising CUTA, IZUMO2, PALB2, and PSCA.

Research has shown that PALB2 and PSCA are closely implicated in DNA repair, radiotherapy response, and tumor progression (28-31). Direct evidence from other reports linking CUTA and IZUMO2 to radiotherapy response remains limited. However, in this study these four genes functioned collectively as an integrated signature to predict radiotherapy response, rather than exerting effects through individual gene functions. This signature was able to effectively predict the OS of patients in both the training and external validation cohorts.

After adjustment for clinicopathological variables, the risk score remained associated with OS, indicating that the association was not entirely explained by factors such as tumor stage or age. Consistent trends in prognostic classification of the model were noted across all the analyzed cohorts; however, additional validation is essential to assess its efficacy in prospective clinical scenarios in RC.

Through functional enrichment analysis, we explored the possible transcriptional regulatory mechanisms relevant to risk stratification. In the high-risk group, gene sets related to immunity and inflammation, as well as pathways related to EMT, showed significant enrichment. Conversely, the low-risk group demonstrated relative enrichment in signatures related to DNA replication and genome maintenance. These results indicated distinct transcriptional profiles among different risk groups. Notably, the co-occurrence of inflammatory and EMT-related pathways in the high-risk group may be associated with cellular plasticity and microenvironmental interactions, while the low-risk group appear to preserve transcriptional characteristics linked to proliferative capacity and genomic stability.

An immunological analysis revealed distinct immune environment characteristics among the different risk groups. High-risk tumors exhibited stronger immune infiltration signals as indicated by ESTIMATE and ImmuneScore analyses, while low-risk tumors showed higher immunogenicity-related scores, including IPS-derived metrics. An analysis of the immune CP genes showed that PDCD1 expression was significantly elevated in the high-risk group, while LAG3 and CTLA4 displayed upward trends. Beyond immune microenvironment differences, somatic mutation profiling revealed distinct genomic landscapes between the risk groups. The high-risk group was characterized by mutations in CACNG3, SOGA1, ZNF521, PCLO, and ROS1, whereas the low-risk group was enriched in PIK3CA and SYNE1 mutations. Notably, no significant differences in the TMB or MSI status were detected between the groups.

Taken together, these integrated immune and genomic analyses suggested that the risk signature was associated with tumor immune microenvironment characteristics and distinct somatic mutation patterns. These multi-dimensional immune and genomic features may contribute to tumor heterogeneity and differential therapeutic responses in RC, providing potential insights for patient stratification and immunotherapy optimization.

A single-cell analysis allowed us to identify the specific cell types in which the prognostic genes were expressed. PSCA expression was largely restricted to epithelial cells, and it marked a subset that was in a specific transcriptional state. Compared with the PSCA− epithelial cells, the PSCA+ epithelial cells were enriched at later pseudotime positions and exhibited broader inferred interactions with immune and stromal populations. At the pathway level, these differences were reflected across multiple ligand-receptor-defined signaling axes identified by CellChat analysis, including representative pathways such as SPP1, MIF, and MDK signaling. These pathways have previously been associated with tumor microenvironment communication. Accordingly, the observed differences are most appropriately interpreted as variations in transcriptional and interaction contexts linked to epithelial heterogeneity, rather than evidence of pathway activation or causal signaling mechanisms.

The study had a number of limitations. All the data were acquired from publicly available clinical databases, and all the analyses were retrospective. The functional characterization of the prognostic module genes was primarily conducted using GSEA, GO, and KEGG analyses, without biological experimental validation. Due to limitations of the available bulk and single-cell transcriptomic datasets, a detailed high-resolution quantification of immune cell infiltration, as well as further analyses of PSCA function in radiotherapy resistance and its co-expression with other signature genes, could not be performed. In addition, while the prognostic signature was derived from radiotherapy-related gene datasets, it was not specifically designed to predict radiotherapy response. Future studies incorporating single-cell sequencing or spatial transcriptomics, together with prospective cohorts and experimental validation, are needed to further elucidate the immune microenvironment heterogeneity and PSCA associated regulatory mechanisms across risk groups.

In summary, this study developed a radiotherapy resistance-related transcriptional signature in LARC that is closely linked to patient prognosis and reflects distinct immune, transcriptional, and epithelial states. Integrative bulk and single-cell analyses suggest that this signature captures key aspects of tumor heterogeneity relevant to radiotherapy response.

Conclusions

We established a novel four-gene prognostic signature for LARC, and the single-cell analysis suggested that PSCA may be a key gene in this signature. Further in vitro and in vivo functional validation is warranted.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0554/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0554/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0554/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. The 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. 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]
2. Airoldi M, Roselló S, Tarazona N, et al. Advances in the management of locally advanced rectal cancer: A shift toward a patient-centred approach to balance outcomes and quality of life. Cancer Treat Rev 2025;140:103015. [Crossref] [PubMed]
3. Kagawa Y, Smith JJ, Fokas E, et al. Future direction of total neoadjuvant therapy for locally advanced rectal cancer. Nat Rev Gastroenterol Hepatol 2024;21:444-55. [Crossref] [PubMed]
4. Bahadoer RR, Dijkstra EA, van Etten B, et al. Short-course radiotherapy followed by chemotherapy before total mesorectal excision (TME) versus preoperative chemoradiotherapy, TME, and optional adjuvant chemotherapy in locally advanced rectal cancer (RAPIDO): a randomised, open-label, phase 3 trial. Lancet Oncol 2021;22:29-42. [Crossref] [PubMed]
5. Conroy T, Bosset JF, Etienne PL, et al. Neoadjuvant chemotherapy with FOLFIRINOX and preoperative chemoradiotherapy for patients with locally advanced rectal cancer (UNICANCER-PRODIGE 23): a multicentre, randomised, open-label, phase 3 trial. Lancet Oncol 2021;22:702-15. [Crossref] [PubMed]
6. Airoldi M, Roselló S, Tarazona N, et al. Advances in the management of locally advanced rectal cancer: A shift toward a patient-centred approach to balance outcomes and quality of life. Cancer Treat Rev 2025;140:103015. [Crossref] [PubMed]
7. Wu Y, Song Y, Wang R, et al. Molecular mechanisms of tumor resistance to radiotherapy. Mol Cancer 2023;22:96. [Crossref] [PubMed]
8. Sharapov MG, Karmanova EE, Gudkov SV. Mechanisms of Cancer Cell Radioresistance: Modern Trends and Research Prospects. Biophysics 2024;69:1064-88.
9. Deng S, Vlatkovic T, Li M, et al. Targeting the DNA Damage Response and DNA Repair Pathways to Enhance Radiosensitivity in Colorectal Cancer. Cancers (Basel) 2022;14:4874. [Crossref] [PubMed]
10. Qin F, Bian Z, Jiang L, et al. A novel high-risk model identified by epithelial-mesenchymal transition predicts prognosis and radioresistance in rectal cancer. Mol Carcinog 2024;63:2119-32. [Crossref] [PubMed]
11. Chen J, Wang S, Ding Y, et al. Radiotherapy-induced alterations in tumor microenvironment: metabolism and immunity. Front Cell Dev Biol 2025;13:1568634. [Crossref] [PubMed]
12. Gai K, Shi X, Xu F, et al. Efficacy and safety of thoracic radiotherapy combined with anti-angiogenic therapy and immunochemotherapy for advanced non-small cell lung cancer patients: a retrospective study. Front Oncol 2025;15:1640306. [Crossref] [PubMed]
13. Beckers C, Pruschy M, Vetrugno I. Tumor hypoxia and radiotherapy: A major driver of resistance even for novel radiotherapy modalities. Semin Cancer Biol 2024;98:19-30. [Crossref] [PubMed]
14. Jin Y, Jiang J, Mao W, et al. Treatment strategies and molecular mechanism of radiotherapy combined with immunotherapy in colorectal cancer. Cancer Lett 2024;591:216858. [Crossref] [PubMed]
15. Grigoriadis D, Tsifintaris M, Giannakakis A, et al. Public Omics Explorer (POE): Enabling integrative semantic search across GEO omics datasets based on PubMed publications. Comput Struct Biotechnol J 2025;27:4802-12. [Crossref] [PubMed]
16. Dong X, Du MRM, Gouil Q, et al. Benchmarking long-read RNA-sequencing analysis tools using in silico mixtures. Nat Methods 2023;20:1810-21. [Crossref] [PubMed]
17. Lin S, Cai K, Chen A, et al. Identification of key modules and hub genes for sepsis-induced myopathy using weighted gene co-expression network analysis. Front Genet 2025;16:1607575. [Crossref] [PubMed]
18. Böge FL, Zacharias HU, Becker SC, et al. Using deep neural networks and LASSO regression to predict miRNA expression changes based on mRNA data. Front Bioinform 2025;5:1566162. [Crossref] [PubMed]
19. Beyene KM, Chen DG. Time-dependent receiver operating characteristic curve estimator for correlated right-censored time-to-event data. Stat Methods Med Res 2024;33:162-81. [Crossref] [PubMed]
20. Luo J, Lu Q, He M, et al. gdGSE: An algorithm to evaluate pathway enrichment by discretizing gene expression values. Comput Struct Biotechnol J 2025;27:1772-83. [Crossref] [PubMed]
21. Zeng D, Ye Z, Shen R, et al. IOBR: Multi-Omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and Signatures. Front Immunol 2021;12:687975. [Crossref] [PubMed]
22. Okada D, Zhu J, Shota K, et al. Systematic evaluation of the isolated effect of tissue environment on the transcriptome using a single-cell RNA-seq atlas dataset. BMC Genomics 2025;26:416. [Crossref] [PubMed]
23. Racle J, Gfeller D. EPIC: A Tool to Estimate the Proportions of Different Cell Types from Bulk Gene Expression Data. Methods Mol Biol 2020;2120:233-48. [Crossref] [PubMed]
24. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
25. Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell 2021;184:3573-3587.e29. [Crossref] [PubMed]
26. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
27. GilisJPerinLMalfaitMDifferential detection workflows for multi-sample single-cell RNA-seq data.bioRxiv 2023;2023.12.17.572043.
28. Foo TK, Xia B. BRCA1-Dependent and Independent Recruitment of PALB2-BRCA2-RAD51 in the DNA Damage Response and Cancer. Cancer Research 2022;82:3191-7.
29. Veenstra CM, Abrahamse P, Hamilton AS, et al. Breast, Colorectal, and Pancreatic Cancer Mortality With Pathogenic Variants in ATM, CHEK2, or PALB2. J Clin Oncol 2025;43:1587-96. [Crossref] [PubMed]
30. Shah Z, Tian L, Li Z, et al. Human anti-PSCA CAR macrophages possess potent antitumor activity against pancreatic cancer. Cell Stem Cell 2024;31:803-817.e6. [Crossref] [PubMed]
31. Stein MN, Dumbrava EE, Teply BA, et al. PSCA-targeted BPX-601 CAR T cells with pharmacological activation by rimiducid in metastatic pancreatic and prostate cancer: a phase 1 dose escalation trial. Nat Commun 2024;15:10743. [Crossref] [PubMed]

(English Language Editor: L. Huleatt)