Differential gene expression analysis was performed using sequencing data from the human colon cancer cells HCT-116 (purchased from the Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai, China), which consisted of three wild-type (HCT116^WT) and three radioresistant (HCT116^RR) cell lines, to identify genes associated with radioresistance. ECM and cell adhesion-related genes were screened, and the survival data from the COAD colon cancer cohort in the TCGA database were used for analysis (Scheme 1).

The TCGA-COAD cohort data were downloaded from the UCSC Xena website (http://xena.ucsc.edu/). The data used in the analysis included gene expression data from 471 tumor tissues, survival data from 454 patients, clinical phenotype data from 478 patients, and tumor mutation data from 399 patients. Risk scores were validated using an external dataset, GSE40967 (7), which was designed to validate whether analytically constructed risk scores could predict prognostic risk in a separate cohort of patients. This dataset contains survival and gene expression assay data for 562 patients.

HCT116 cells were exposed to a linear accelerator (Varian Clinic 21 ES; Varian Medical Systems, Crawley, UK), which generates high-energy X-rays (6 MeV) at 0.99 Gy/min. Radioresistant cell lines (HCT116 RR) were developed from the parental cell line (HCT116) by increasing the X-ray dose of fraction irradiation stepwise from 0.5 to 2 Gy/day in vitro (2). Cells were initially exposed to 0.5 Gy/day of X-rays for 5 days. Subsequently, the cells were exposed to 1 Gy/day of X-rays for 10 days. Thereafter, the surviving cells were exposed to 1.5 Gy/day of X-rays for 15 days. Cells that could proliferate under exposure to 1.5 Gy/day of X-rays were further exposed to 2 Gy/day of X-rays. If these cells proliferated constantly under exposure to 2 Gy/day of X-rays for more than 30 days, it was determined that a HCT116 RR had been obtained.

Total RNA was extracted from the HCT116 cell line using TRIzol^® Reagent according the manufacturer’s instructions (Magen). RNA samples were detected based on the A260/A280 absorbance ratio with a Nanodrop ND-2000 system (Thermo Scientific, USA), and the RIN of RNA was determined by an Agilent Bioanalyzer 4150 system (Agilent Technologies, CA, USA). Only qualified samples were used for library construction.

Paired-end libraries were prepared using a ABclonal mRNA-seq Lib Prep Kit (ABclonal, China) following the manufacturer’s instructions. The mRNA was purified from 1 μg of total RNA using oligo (dT) magnetic beads followed by fragmentation carried out using divalent cations at elevated temperatures in ABclonal First Strand Synthesis Reaction Buffer. Subsequently, first-strand cDNAs were synthesized with random hexamer primers and Reverse Transcriptase (RNase H) using mRNA fragments as templates, followed by second-strand cDNA synthesis using DNA polymerase I, RNAseH, buffer, and dNTPs. The synthesized double stranded cDNA fragments were then adapterligated for preparation of the paired-end library. Adaptor-ligated cDNA was used for PCR amplification. PCR products were purified (AMPure XP system), and library quality was assessed on an Agilent Bioanalyzer 4150 system. Finally, the library preparations were sequenced on an Illumina Novaseq 6000 (or MGISEQ-T7), and 150-bp paired-end reads were generated.

Raw data in fastq format were firstly processed through in-house perl scripts. In this step, the adapter sequence was removed and low quality reads were filtered out (the number of lines with a string quality value less than or equal to 25 accounted for more than 60% of the entire reading and the N ratio (i.e., base information could not be determined) was greater than 5% of total reads) to obtain clean reads that could be used for subsequent analysis. Then, the clean reads were separately aligned to the reference genome in the orientation mode using the HISAT2 software (http://daehwankimlab.github.io/hisat2/) to obtain mapped reads. FeatureCounts (http://subread.sourceforge.net/) was used to count the read numbers mapped to each gene. Then, the FPKM of each gene was calculated based on the length of a gene and the read count mapped to that gene. Differentially expressed genes (DEGs) were identified using the “DESeq2” R package. Adjusted p-value < 0.05 and |log₂FC| > 1 were used as cut-off points to identify DEGs for subsequent analysis.

The candidate genes were used to generate prognosis-related risk scores. The patients in the TCGA-COAD cohort with both gene expression information from tumor samples and complete survival information served as the training set. Univariate Cox regression analysis was used to screen for DEGs in the training set that were strongly associated with patient survival, using a p-value < 0.05. Next, LASSO regression was applied to further establish radioresistance-related risk profiles. The prognostic risk score model was established with the following formula: risk score = expression level of Gene1 × β1 + expression level of Gene2 × β2 +…+ expression level of Genen × βn (where β is the regression coefficient calculated by the LASSO regression). Risk scores were calculated for each patient using a risk-score model. The samples were assigned to high-risk or low-risk groups according to the median risk score. Kaplan–Meier curves were used to compare the differences in overall survival (OS) between the high-risk and low-risk groups. ROC curves for the 1-, 3-, 5- and 7-year OS were generated for the two groups. The established risk score was evaluated in an external independent set, GSE40967, to assess its performance in prognosis prediction.

Independent prognostic factors selected by univariate and multivariate Cox regression analyses were used to construct nomograms to assess the probability of OS. Subsequently, calibration curves were used to estimate whether the predicted survival outcome was close to the actual outcome. DCA was performed to confirm the clinical utility of the nomogram model.

Analyses were performed in the R software (R software, version 4.1.0). The R and package versions can be found in Supplementary Figure S2 . The Wilcoxon test was used for statistical comparisons between two groups, the Kruskal–Wallis test for statistical comparisons between multiple groups. The Bonferroni method was used to correct for multiple comparisons. An adj.p-value < 0.05 was regarded as statistically significant. Cox regression analysis was performed for outcomes adjusted with covariates for which the p-value was <0.05 in a univariable Cox analysis. the Spearman correlation coefficient for correlations between two continuous variables, except for group comparisons specific to high-throughput data.

We first verified the radioresistance of the HCT116^WT and the HCT116^RR cell line using clone formation assay ( Supplementary Figure S1 ). We first screened for DEGs using the HCT-116 radiotherapy-tolerant cell line HCT116^RR and the HCT116^WT cell line as comparative samples. A total of 493 DEGs (RR vs. WT) were identified, of which 259 genes were upregulated and 234 genes were downregulated ( Supplementary Table S1 ), using a p < 0.05 and |log₂FC|>1 as thresholds ( Figure 1A ). The top-15 upregulated and downregulated genes were arranged according to the p-value and visualized using a heat map ( Figure 1B ).

GO enrichment analysis ( Figures 1C–E and Supplementary Table S2 ) showed that the DEGs were mainly enriched in the biological processes “extracellular matrix organization”, “extracellular structure organization”, and “external encapsulating structure organization”. In cellular the components and molecular function categories, the DEGs were mainly enriched in the functions “collagen-containing extracellular matrix”, “focal adhesion”, “cell-substrate junction”, “extracellular matrix structural constituent”, and “extracellular matrix binding”, which indicates that the differential genes contribute the most to the functions of extracellular matrix structural constituents and extracellular matrix binding.

Similarly, KEGG enrichment results ( Figure 1F and Supplementary Table S3 ) showed that the DEGs contained the two pathways of focal adhesion and ECM-receptor interaction. Combining the GO and KEGG enrichment results, we determined that the DEGs were mainly enriched for cell adhesion and ECM-related functions and pathways. Therefore, the cell adhesion and ECM-related DEGs were screened as candidate genes for subsequent analyses.

We cross-referenced the cell adhesion and extracellular matrix-related genes with the above-mentioned DEGs by releasing all genes containing the following annotations or their subdivisions: GO:0031012 extracellular matrix; GO:0030198 GO:0031012 extracellular matrix; GO:0030198 extracellular matrix organization; GO:1903053 regulation of extracellular matrix organization; GO:0035426 extracellular matrix-cell signaling; GO:0007155 cell adhesion; GO:0030155 regulation of cell adhesion; GO:0007160 cell-matrix adhesion; and GO:0050840 extracellular matrix binding. A total of 94 overlapping genes from DEGs and ECM-related genes were retained as candidate genes for subsequent analysis ( Figure 1G ). For these 94 genes, the interaction network of the encoded proteins was mapped using the STRING database ( Figure 1H ).

The relationship between candidate gene expression and the OS of patients was analyzed using a one-way Cox regression model. Five genes, TNFRSF13C, CD36, ANGPTL4, LAMB3, and SERPINA1, were found to be associated with OS, as shown in Figure 2A . LASSO regression analysis was performed using these five prognostic genes ( Figure 2B ), and the best model parameter was λ=0.001838643. The linear combination of the five genes and their coefficients in the model were used as the risk score, and the coefficients of the five genes are shown in Table 1 .

Risk scores were calculated for patients in the TCGA COAD cohort, and patients were divided into two groups based on the median risk score: high and low risk ( Figure 2C ). The risk score distribution and associated survival status implied that there were more dead patients with an increase in the risk score. Furthermore, survival analysis showed that there was a significant difference between the low- and high-risk groups (p-value < 0.0001, Figure 2D ). After statistically testing whether there were differences in risk scores among patients with different tumor stages and different sexes, we found that patient risk scores increased significantly with the stage of TNM, the key indicator of tumor progression ( Figure 2E ). However, patient risk scores were independent of sex.

Then, the performance of the risk scores for predicting the prognostic risk of patients was validated using an external independent validation set. We divided patients into high-risk and low-risk groups based on the median risk score as the cut-off point ( Figure 3A ). The Kaplan-Meier survival curve results showed that patients with lower risk scores exhibited higher survival ( Figure 3B ), and there was a significant difference between the high- and low-risk groups (p=0.0054).

To further investigate the differences in functions and pathways between the tumors of the high- and low-risk groups, we performed gene set enrichment analysis on the data. The pathways represented by different lines of red are molecular pathways or functions that were significantly activated in the tumor tissues of patients in the high-risk group, while the pathways represented by different lines of blue are molecular pathways or functions that were significantly activated in the tumor tissues of patients in the low-risk group.

Regarding biological processes, samples from the high-risk group were enriched in “Collagen fibril organization,” “antigen processing and presentation of peptide or polysaccharide antigen via MHC class II,” and “B cell receptor signaling pathway.” The samples from the low-risk group were enriched in “mitochondrial gene expression” and “mitochondrial respiratory chain complex assembly” ( Figure 4A ). Cellular Component analysis showed that “collagen trimer”, “protein complex involved in cell adhesion”, and “collagen containing extracellular matrix” were enriched in high-risk samples, whereas “large ribosomal subunit”, “small ribosomal subunit”, and “organellar ribosome” were enriched in low-risk samples ( Figure 4B ).

Regarding molecular functions, the high-risk group mainly showed enrichment in “collagen binding” and “extracellular matrix structural constituent”, while “oxidoreductase activity acting on NADPH”, “RNA binding”, and “oxidoreduction driven active transmembrane transporter activity” were enriched in the low-risk group ( Figure 4C ). Meanwhile, the KEGG gene sets showed enrichment in “ECM receptor interaction”, “focal adhesion”, and “cell adhesion molecules cams” in the high-risk group, while “ribosome”, “aminoacyl tRNA biosynthesis”, and “peroxisome” were enriched in the low-risk group ( Figure 4D ).

The R package “ESTIMATE” was used to analyze the correlation between immune/stroma scores and risk scores. The results showed that both the immune score ( Figure 5A ) and stroma score ( Figure 5B ) were significantly higher in the high-risk group than in the low-risk group (both p<0.05). In addition, the ssGSEA method was used to estimate the infiltration of 28 immune cell species in samples from the high- and low-risk groups. The infiltration of 18 immune cell species significantly differed between the high- and low-risk groups ( Figure 5C ).

We also compared the expression of 34 immune checkpoint genes and 23 HLA family genes between the high- and low-risk groups and found that the expression of 22 immune checkpoints and 14 HLA family genes significantly differed between the two groups ( Figures 6A, B ). Among them, the top-three immune checkpoint genes and HLA family genes with the most significant differences were CD134, FOXP3, GEM, B7H5, CD134L and HLA-DQAZ, HLA-DMA, HLA-F, HLA-DMB, HLA-DPA1, respectively. To determine the potency of the risk score as a biomarker for predicting drug response in radiotherapy-tolerant patients with CRC, we assessed the drug-sensitivity of patients in the high- and low-risk groups to different antitumor drugs. The results show that a total of 44 drugs showed significant differences in sensitivity between the two groups. Among them, The ten most sensitive drugs in the high-risk group were Bortezomib, Epothilone B, Elesclomol, BEZ235, CGP-60474, Gemcitabine, QL-VIII-58, AZD7762, rTRAIL, and MG-132 ( Supplementary Table S4 ). In terms of response to immunotherapy, patients in the low-risk group had lower TIDE dysfunction scores ( Figure 6B ); however, TIDE exclusion scores did not significantly differ from those in the high-risk group ( Figure 6D ).

To investigate the mechanisms associated with poor prognosis in radiotherapy-tolerant CRC, we analyzed somatic mutations of all genes in the TCGA database. As shown in the waterfall figure ( Figure 7A ), APC, TP53, and TTN were the three most mutated genes in the COAD cohort. KRAS and MUC16 were also identified as mutated genes in high-risk samples ( Figure 7B ). FAT4 and OBSN were representative mutated genes in low-risk samples ( Figure 7C ). When comparing the mutation frequencies between the low- and high-risk group samples, more somatic mutations were observed in the high-risk group, but there was no significant difference (p=0.92) between the two groups ( Figure 8A ). The mutation frequencies of the 26 genes differed between the high- and low-risk groups ( Figure 8B ). Additionally, we visualized the mutation distribution of the five model genes; among them, LAMB3 was the most mutated gene in the sample, and missense mutations were the most common mutation classification ( Figure 7D ).

Risk score, age, and TNM stage were identified as significant predictors of prognosis, using univariate and multivariate Cox regression analyses ( Figure 9 ). To predict the prognosis of each patient, a forest plot was created by integrating the risk score, age, and TNM stage, and we found that the hazard ratio of the risk score was higher than that of the T and N stages ( Figures 9A, B ). The calibration curve showed a high agreement between the risk score prediction and actual survival probability ( Figures 9C–E ). Next, we plotted ROC curves to examine the specificity and sensitivity of this prognostic model ( Figures 9F–H ). The time-dependent ROC curve showed that the AUC of the resulting prognostic model (based on the five-gene risk score and other clinical parameters, including age, sex, and tumor stages) was 0.79, 0.77, and 0.78 at 1, 2, and 3 years, respectively. The calibration curve showed high agreement between the risk score prediction and actual survival probability. In addition, the net benefit of this prognostic model was high, as shown by the DCA curve ( Figures 9I–K ).