Transcriptome sequencing data and clinical information on CC patients were downloaded from The Cancer Genome Atlas (TCGA) (https://www.cancer.gov/tcga/). Patients with incomplete clinical pathological information and missing survival data were excluded. We finally followed up 436 CC patients. We downloaded RNA-seq TCGA transcriptome data in a standardized fragment format per kilobase exon model per million mapped reads. Supplementary Table S1 shows the clinical information of all patients.

To obtain the lncRNA expression profiles, we mapped TCGA sequencing data to lncRNA annotation files in the GENCODE database (http://www.gencodegenes.org) and eliminated lncRNAs with no expression in > 20% of samples. We screened PCD-related lncRNAs in TCGA-Colon Adenocarcinoma (COAD) data using Spearman correlation analyses based on correlation coefficients > 0.3 and P < 0.05 ( Supplementary Table S2 ). The top 10 lncRNAs for differential PCD genes were visualized using Cytoscape to explore the prognostic value of PCD-related lncRNAs, and those with prognostic value were screened using univariate Cox regression analysis. The expression of PCD-related lncRNAs in CC and normal tissues was analyzed using Wilcoxon signed-rank tests.

We used the ConsensusClusterPlus package in R to determine the number of consistent clusters based on the expression PCD-related prognostic lncRNAs and their potential biological characteristics. Prognostic value and clinicopathological characteristics in subgroups were analyzed using Kaplan–Meier (KM) survival analysis and log-rank tests.

Differential genes among subgroups were screened using the limma package in R with the following criteria: |logFC| ≥ 0.5 and a false discovery rate (FDR) < 0.05 ( Supplementary Table S3 ). Differentially expressed genes (DEGs) among subgroups were analyzed using Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG), clusterProfiler, org.Hs.eg.db, and the enrichplot package in R to determine their biological functions. Differences in tumor signal pathways between patients in clusters 1 and 2 were screened via Gene Set Enrichment Analysis (GSEA) using FDR < 0.05 as the criterion.

We calculated immune and stromal scores in CC to compare immune infiltration between the subgroups using the ESTIMATE algorithm. The abundance of immune cells was determined using TIMER (19), CIBERSORT (20, 21), QUANTISEQ (22), Microenvironment Cell Populations-counter (MCPCOUNTER) (23), XCELL (24), and the Estimating the Proportion of Immune and Cancer cells (EPIC) algorithm (25). Immune pathways in groups were quantified using the GSVA package in R. Differences in immune checkpoints between subgroups were evaluated using Wilcoxon tests.

We constructed a prognostic model by screening eight PCD-related lncRNAs and then optimized readability using the LASSO regression algorithm based on their expression and corresponding risk coefficients. The formula was calculated as follows:

The training and validation sets were divided into high- and low-risk groups based on the median risk score. The overall survival (OS) of patients was assessed using KM curves and log-rank tests. The effectiveness of the model was evaluated using receiver operator characteristic (ROC) curves. Whether risk scores and clinical characteristics were independent prognostic factors for patients were evaluated using univariate and multivariate Cox analyses. Differences in risk scores, immune scores, and clinicopathological characteristics were assessed using Wilcoxon tests. We used the RMS package in R to integrate clinical characteristics with risk scores and construct a nomogram and calibration curve to realize a quantitative prognostic tool. Consistency between the prognostic assessment by the nomogram and actual results was evaluated using a calculation with a guide method that included 1,000 resamples. The diagnostic power of the nomogram and individual predictors was evaluated using ROC curves. The clinical benefits conferred by prognostic evaluation of the nomogram and a single predictor were further compared using decision curve analysis (DCA).

Somatic CC mutation data were retrieved from TCGA (https://www.cancer.gov/tcga/). The TMB was calculated by dividing the total number of mutations by the size of the coding region of the target gene. Patients were classified as having high or low TMB based on the median value. We also visualized the top 20 genes with the highest mutation frequency in high- and low-risk groups using the maftools package in R. Copy number variations (CNV) between subgroups were calculated using Chi-square tests, and the positions of CNVs on the chromosome were visualized using the Rciorcos package in R.

The sensitivity of subgroups of patients to immunotherapy was assessed using http://tide.dfci.harvard.edu. Based on the Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org), we used the pRRophetic package in R to evaluate the sensitivity of patients to chemotherapeutic agents and visualized the three-dimensional (3D) molecular structure of each agent using the cMAP database (https://clue.io/cmap).

We collected CC tumor and adjacent cancerous tissue from patients who underwent CC resection between 2016 and 2021 at the Colorectal Cancer Center of Shanghai Tenth People’s Hospital. After excluding patients with incomplete clinical data, 96 were included in the follow-up analysis. The study was approved by the Ethics Committee at Shanghai 10th People’s Hospital.

Total RNA from tissues or cells was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and quantified using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). Total RNA was reverse transcribed to cDNA using the PrimeScript™ RT Reagent Kit (Takara Bio Inc., Kusatsu, Japan). Amplicons were quantified using a 7500 Fast Real-time PCR system (Applied Biosystems, Piscataway, NJ, USA). The relative expression (fold change) of the target genes was determined using the 2-ΔΔCT method. Glyceraldehyde 3-phosphate dehydrogenase was the internal control. Supplementary Table 4 shows the sequences of the RT-qPCR primers (Beijing Qingke Biotech Ltd., Beijing, China).

Human HCT-116 and SW480 CC cell lines (American Type Culture Collection, Manassas, VA, USA) were cultured in Dulbecco’s modified Eagle’s medium containing 10% fetal bovine serum (FBS; both from Thermo Fisher Scientific Inc., Waltham, MA, USA) with 1% streptomycin and penicillin, at 37°C under a humidified 5% CO₂ atmosphere.

Cells were transiently transfected using Lipofectamine 2000 (Invitrogen) as described by the manufacturer. Small interfering RNA (siRNA) was purchased from GenePharma (Shanghai, China) The sense and anti-sense si-U62317.4 sequences were: 5′-GAAGAGAAGGACAAGUUGACG-3′ and 5′-UCAACUUGUCCUUCUCUUCUG-3′, respectively. We created a stable knockout U62317.4 cell line using short hairpin RNA (shRNA) targeting U62317.4 (sh-U62317.4; General Biosystems, Anhui, China) with the sense and anti-sense sequences, 5′-GATACTTGATCCTGATAAA-3′ and 5′- TTTATCAGGATCAAGTATC-3′, respectively.

Lentiviral particles were obtained as described (16). HCT-116 cells were directly infected with Polybrene (Santa Cruz Biotechnology, Dallas, TX, USA) for 24 h and then, transfected cells were screened for 7–10 days with 2 μg/mL of puromycin (Invitrogen).

We transfected HCT-116 or SW480 cells for 48 h and then seeded 1 × 10³ cells/well in 96-well plates. After 24, 48, and 72 h, the cells were incubated with CCK-8 reagent (10 µL) in serum-free medium (100 µL) for 2 h at 37°C. The optical density at 450 nm was measured using a microplate reader. Similarly transfected CC cells were seeded in six-well plates and cultured in DMEM medium containing 10% FBS for 7 days to analyze colony formation. The cells were then fixed with paraformaldehyde and stained with 0.1% crystal violet. Colonies were counted using ImageJ software.

Cell invasion and migration were measured using Transwell chambers (Corning CoStar, Tewksbury, MA, USA). Cells (4 × 10⁴) were resuspended in serum-free medium and added to the upper part of Matrigel^®-coated chambers. Medium (500 µL) containing 10% FBS was added to the lower chamber. The cells were fixed with paraformaldehyde 48 h later and stained with 0.5% crystal violet. Cells that passed through the bottom of the membrane were assessed and counted using a microscope. The cells were resuspended in serum-free medium and seeded in six-well plates until they reached 100% confluence, when they were damaged by being gently scratched with a 200-μL sterile micropipette tip for wound healing assays. Representative images of cell migration were acquired after 0, 24, and 48 h, and the cell migration rate was determined by time-lapse analysis using ImageJ software.

We assayed apoptosis using Flow Apoptosis Kits (BD Biosciences, San Diego, CA, USA) as described by the manufacturer. Briefly, digested cells were washed twice with cooled PBS and stained using FITC Annexin V apoptosis detection kits (BD Biosciences); thereafter, cell populations were evaluated using a BD LSRFortessa™ analyzer (BD Biosciences).

Four-week-old female BALB/C nude mice were randomly assigned to two groups (n = 5 per group) to determine tumor formation in vivo. CC cells stably transfected with sh-normal control (NC) or sh-U62317.4 were implanted subcutaneously into the axillae of nude mice. One week after injection, tumor volumes (cm³) were measured every 3 days as (length × width²)/2. The mice were euthanized 21 days later, and tumors were removed and weighed. The Animal Experiment Ethics Committee of Tongji University approved all the animal experiments.

Proteins were collected from CC cells in RIPA lysis buffer (Invitrogen) containing phenylmethylsulfonyl fluoride (PMSF; Bio-Rad Laboratories Inc., Hercules, CA, USA). Proteins were resolved by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and then transferred onto polyvinylidene difluoride membranes (Invitrogen). Non-specific antigen binding was blocked by shaking the membranes in 5% skim milk for 2 h at 37°C. Thereafter, the membranes were incubated with primary antibody overnight at 4°C, followed by secondary antibody for 2 h. The blots were visualized and analyzed using the Odyssey system (LI-COR Biosciences, Lincoln, NE, USA). Supplementary Table 5 shows details of the antibodies.

Data were statistically analyzed, and images were generated using R and GraphPad Prism 8.03 (GraphPad Software Inc., San Diego, CA, USA). The Wilcoxon signed-rank and Kruskal-Wallis tests were used for between-group comparisons. Values with P < 0.05 were considered statistically significant.

Figure 1 shows a flow chart of the study. We initially re-annotated TCGA expression matrix to distinguish between mRNAs and lncRNAs and identified 199 genes related to PCD among which, 33, 52, and 114 were associated with apoptosis (26), pyroptosis (9, 27), and ferroptosis (28). Figures 2A–C shows that the expressions of most PCD-related genes significantly differed between CC and normal tissues. We then screened the lncRNAs that were the most closely associated with the expression of these genes using correlation analysis. Among the 437 PCD-related lncRNAs, 85, 161, and 297 were associated with pyroptosis, apoptosis, and ferroptosis, respectively ( Figures 2D–F ). Univariate Cox analysis of the screened PCD-lncRNAs identified 41 PCD prognosis–related lncRNAs ( Figures 3A, B ). Differential analysis showed that the expression of these lncRNAs significantly differed between tumor and normal tissues ( Figure 3C ). Consensus clustering based on the expression of these lncRNAs revealed biological differences among different CC subgroups. Combining the consensus matrix cumulative distribution function curve and the delta area plot showed that interference between clusters was minimal and the classification was significant when K = 2 ( Figures 3D–F ). KM survival curves showed that the OS rate was significantly better for cluster 1 than that cluster 2 ( Figure 3G ). We then constructed a heat map to compare the clinical characteristics between these subgroups and found high scores for T stage, lymphatic metastasis, and tumor stage in cluster 1 ( Figure 3H ).

The survival results revealed evident prognostic differences between the CC subgroups. To explore the potential reason for this, we analyzed DEGs in the subgroups using GSEA, GO, and KEGG. The GSEA results showed that the immune-related B cell receptor, chemokine, JAK-STAT, T cell receptor, and Toll-like receptor signaling pathways were abnormally activated in cluster 1 ( Figure 4A ). The GO results showed that the DEGs were involved in several immune-related functions, which was consistent with the GSEA findings ( Figures 4B–D ). The KEGG analysis showed that the DEGs were primarily involved in cell growth, chemokines, cytokines, and phagosomes related to signaling pathways ( Figure 4E ).

Considering that the DEGs were primarily related to immune-related functions, we examined the immune cell distribution in the immune microenvironment of CC. The amounts of immune cell infiltration differed between clusters 1 and 2 ( Figure 5A ). Cluster 1 had abundant memory CD4 T cells and eosinophils, whereas cluster 2 had more follicular helper T and resting natural killer (NK) cells ( Figures 5B–E ). We further explored differences in the distribution of immune and stromal ratios between the two clusters using the ESTIMATE algorithm. Figures 5F–H shows that cluster 2 had lower immune, stromal, and ESTIMATE scores than those of cluster 1. A comparative analysis revealed fully suppressed immune function in cluster 2 ( Figure 5I ). Considering the importance of immune checkpoints in immunotherapy, we analyzed the expression of immune checkpoints in the two clusters and found that all immune checkpoint genes except the TNF Receptor Superfamily Members (TNFRSF) 25 and 14 were significantly suppressed in cluster 2 ( Figure 5J ). Thus, the lower OS rate of patients in this cluster might have been due to suppressed immune function.

LASSO Cox analysis of PCD prognosis–related lncRNAs in TCGA identified the filter genes AC004080.1, AC078923.1, AC114730.3, AC156455.1, long intergenic non-protein coding RNA 1419 (LINC01419), LINC00997, NF-kappaB interacting lncRNA (NKILA), and U62317.4 with which to construct a prognostic model. Risk scores were calculated based on the regression coefficients and expression of these eight genes ( Figures 6A–C ). The risk formula was 0.5023 * expression (U62317.4) + 0.3836 * expression (AC078923.1) + 0.0082 * expression (AC114730.3) + 0.075 * expression (LINC00997) + 0.1341 * expression (NKILA) + 0.1278 * expression (AC156455.1) + 0.0195 * expression (LINC01419) + 0.1131 * expression (AC004080.1). To improve the accuracy and effectiveness of the prognostic model, we randomly assigned the patients to training (N = 260) and test (N = 176) groups and then divided them into high-and low-risk groups according to the median risk score. The KM survival curves indicated worse OS for patients in the training group with high risk scores than those with low risk scores ( Figure 6D ). The risk curves indicated a significantly higher death rate among patients with high risk scores than those with low risk scores ( Figure 6G ). The ROC curves showed that the prognostic model accurately predicted the 1, 3, and 5-year survival rates of CC patients ( Figure 6J ). This model was further verified in the test group ( Figures 6E, H, K ) and the entire TCGA database ( Figures 6F, I, L ). These findings indicated that our prognostic model is unbiased and can be used as a reference tool for predicting the OS rate of CC.

We further confirmed the independence of the CC prognostic model using univariate and multivariate Cox analyses. Figures 7A, B shows that age (hazard ratio [HR], 1.042; 95% confidence interval [CI], 1.019–1.067; P < 0.001) and risk score (HR, 1.493; 95% CI, 1.346–1.655; P < 0.001) were independent prognostic factors for CC patients. We then evaluated the relationship between the risk score and PCD prognosis–related lncRNA. A heatmap showed significantly increased expression of all eight PCD-related lncRNAs in the high-risk group ( Figure 7C ). Risk scores were statistically different in cluster stratification, tumor stage, and lymph node metastasis ( Supplementary Figure 1A ). In addition, risk scores could assess the prognosis of patients in multiple clinical subgroups, except for patients in the T1-2 group ( Supplementary 1B ). We constructed a nomogram that included patient age, sex, tumor stage, distant metastasis, lymph node metastasis, and risk score to accurately quantify survival rates ( Figure 7D ). The calibration curves showed that the actual 1-, 3-, and 5-year OS rates of patients and those estimated by the nomogram were close ( Figure 7E ). The areas under the ROC curve (AUC; Figure 7F ) revealed that the 1, 3, and 5-year survival rates determined by the nomogram were accurate (AUC = 0.820 0.824, and 0.838, respectively). The DCA decision curve showed that the net rate of return for the OS rates assessed by the combined model was better than other clinical characteristics ( Figures 7G–I ). The ROC curve also showed that the combined model was more sensitive than other clinical features ( Figures 7J–L ). These results showed that our nomogram can help clinicians plan accurate follow-up strategies.

We explored the potential biological functions of and signaling pathways enriched by the DEGs in the high- and low-risk groups ( Supplementry Table 6) . The GO results showed that DEGs were significantly enriched in the processes of epithelial-mesenchymal transition, Wnt signaling pathway regulation, and cytokine stimulation ( Supplementary Figure 2A ). The KEGG results revealed that DEGs were primarily involved in Wnt, Hippo, extracellular matrix-receptor interactions, transforming growth factor-beta, and cytokine-cytokine receptor interaction signaling pathways ( Supplementary Figure 2B ).

Figure 8A shows an immune infiltration heat map based on the TIMNER, CIBERSORT, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC algorithms. Among immune functions, difference analysis showed that Apcin (APC) co-inhibition, APC co-stimulation, checkpoint, cytolytic activity, inflammation promotion, MHC class I, T cell co-inhibition, and T cell co-stimulation significantly differed between the high-and low-risk groups ( Figure 8B ). These findings showed that the eight PCD-related lncRNA prognostic characteristics of CC are somewhat related to immune cell infiltration.

In view of the evident prognostic differences in the high- and low-risk patients, we analyzed the TMBs and CNVs to further uncover underlying causes. Figure 8C shows lower TMB levels in patients with high risk scores than those with low risk scores. Further correlation analysis showed that the risk score negatively correlated with TMB ( Figure 8D ). We then divided the patients into groups with a high or low TMB based on the median TMB. KM curves revealed a significantly worse OS rate for patients with a high TMB than those with a low TMB ( Figure 8E ). Considering that the risk score and the TMB have good prognostic value in CC patients, we evaluated the synergistic effects of these scores on the prognostic stratification of CC. We found that the TMB did not affect the assessment of the prognosis of CC patients according to risk score. Survival in risk score subtypes significantly differed between groups with a high or low TMB ( Figure 8F ). These results indicated that the risk score might be an effective indicator that can evaluate prognosis independently of the TMB. We assessed the distribution of somatic variations in CC driver genes between the high- and low-risk groups. A waterfall chart shows the top 20 genes with the highest mutation frequency ( Figures 8G, H ). Our analysis of mutation annotation files in TCGA cohort revealed that FAT4, OBSCN, PCLO, ABCA13, ZFHX4, DNAH11, RYR2, and USH2A significantly differed between the high- and low-risk groups ( Supplementary Table 7 ). In addition, we performed CNV analysis of the top 40 genes with the greatest differences between the high- and low-risk groups. We found a higher frequency of CNV mutations in DEGs between the high-and low-risk groups, where CNV expansion occurred in SALL4, KRT16, PEG10, EPHB6, PRR9, GP2, FGF19, and SULT1C4, whereas CST1, H2BC8, HOXD10, HOXD10, HOXD11, CLCA4, DKK1, IGFL1, ALDOB, CALCA, and PLIN4 underwent CNV deletion ( Figure 8I ). The locations of CNV mutations in these DEGs are shown in Figure 8J . These results may provide new ideas for studying gene mutations in PCD-related lncRNAs in CC.

We evaluated the efficacy of immunotherapy in the high- and low-risk groups using the website http://tide.dfci.harvard.edu/. A higher tumor immune dysfunction and exclusion (TIDE) score implies a higher possibility of immune escape, indicating that the patient has a lower benefit from immunotherapy. As shown in Figure 8K and Supplememtary Table 8 , the TIDE scores of patients in the high-risk group were significantly higher than those in the low-risk group, implying that the low-risk group patients would benefit from immunotherapy. Patients in the high-risk group had lower microsatellite instability (MSI) scores ( Figure 8L ), whereas patients in the low-risk group had higher T cell exclusion scores ( Figure 8M ). T cell dysfunction did not significantly differ between the subgroups ( Figure 8N ). These results support a basis for novel customized immunotherapies for CC patients.

Considering that chemotherapy and targeted therapies are popular strategies for treating CC, understanding the sensitivity of subgroups of patients to such drugs is important. We predicted the sensitivity of the high-and low-risk groups to agents that are commonly administered to CC patients. Figures 9A–F shows that the group with low-risk scores was more sensitive to PD.0325901 METFORMIN MK.2206 AZD8055 PD.0332991, and sorafenib, whereas that with high-risk scores was more sensitive to imatinib, lapatinib, PHA.665752, and MS.275 ( Figures 9G–J ). We determined the 3D structure of four chemotherapy drugs that could be used for patients in the high-risk group using the CMAP database ( Figures 9K-N ). These results should facilitate the application of precise and or new personalized medicines for treating CC.

We selected 96 CC patients from the Colorectal Cancer Center of Shanghai Tenth People’s Hospital to verify the external dataset and examined the expression of the eight PCD-associated lncRNAs in CC and adjacent tissues by RT-qPCR. The expression trends were consistent with those listed in TCGA ( Figure 10A ). We calculated the risk scores of the patients and divided them into high- and low-risk groups according to the median risk score. KM curves showed that the OS was significantly worse in the high- risk group than that in the low-risk group ( Figure 10B ). The AUC showed that the risk model could predict the outcome of CC patients ( Figure 10C ). Univariate and multivariate Cox regression analyses of clinical characteristics and risk scores determined that risk scores, stage and age were independent prognostic factors ( Figures 10D, E ).

The findings that the lncRNA U62317.4 had the highest risk coefficient in the prognostic model indicated that U62317.4 is closely associated with the prognosis of CC. However, its role in the occurrence and development of CC is unclear. Therefore, we explored whether U62317.4 is involved in the malignant progression of CC. We synthesized a siRNA against U62317.4 and evaluated its interference via RT-qPCR. The expression of U62317.4 in HCT-116 and SW480 cells was significantly reduced after transfection with si-U62317.4 ( Figures 11A, B ). Silencing U62317.4 significantly inhibited CC cell viability and clone formation ( Figures 11C–E ) and significantly increased the apoptosis of CC cells ( Figures 11F, G ). Metastasis is another important feature of malignant tumors. Thus, we explored the role of U62317.4 in CC metastasis using wound healing and Transwell assays. Knockdown of U62317.4 significantly reduced wound healing ability as well as the migratory and invasive capacity of CC cells ( Figures 11H, J ). We constructed a nude mouse subcutaneous tumor model. to evaluate the tumorigenicity of U62317.4 in vivo. Silencing U62317.4 significantly inhibited the growth of CC in these mice ( Figures 11K–M ). Western blotting revealed changes in the abundance of proliferation- and apoptosis-related proteins when U62317.4 was silenced. Figures 11I shows that the expression of ki-67, PCNA, and BCL-2 decreased, whereas that of cleaved caspases-3 (cleaved casp-3) and -8 increased when U62317.4 was silenced in CC. These results showed that silencing U62317.4 inhibits CC cell proliferation and metastasis.