A gene signature related to programmed cell death to predict immunotherapy response and prognosis in colon adenocarcinoma

Data collection

Single nucleotide variation (SNV) mutation information, the expression data in transcripts per million (TPM) format, and clinical follow-up information of COAD patients were collected from the TCGA database (https://portal.gdc.cancer.gov/) (Shahrajabian & Sun, 2023). The RNA-Seq expression profile was log2-transformed and Ensembl was converted to gene symbol. The SNV mutation information was processed by mutect2. A sum of 427 COAD samples with survival time longer than 30 days were obtained (Table S1). In addition, the cBioportal database (https://www.cbioportal.org/) was accessed to acquire the RNA-seq expression and survival data of COAD in AC-ICAM cohort (validation set), and a total of 341 tumor samples with a survival time longer than 30 days were retained. Based on a previous study, 1,254 PCD-associated genes were extracted for our follow-up study (Zou et al., 2022).

Consensus clustering analysis

Firstly, prognostically relevant genes were selected from a total of 1,254 PCD-related genes by univariate Cox analysis. Then, consensus clustering was conducted using ConsensusClusterPlus R package (parameters: clusterAlg = “pam” and distance = “pearson”) (Chen et al., 2022). Sampling was repeated 500 times at a sampling ratio of 80% each time. The threshold was between two and 10 to decide the optimal number for the clustering.

Evaluation of immune microenvironment and immunotherapy response

CIBERSORT (https://cibersort.stanford.edu/) was used to quantify the density of 22 types of immune cells in COAD samples (Chen et al., 2018). MCP-counter was used to quantify the abundance of three types of stromal cells and eight types of immune cells based on gene expression matrix (Becht et al., 2016).

Immunotherapy benefits to different risk groups were predicted applying TIDE (http://tide.dfci.harvard.edu/) and aneuploidy scores. Potential immune checkpoint blockade therapy responses in COAD patients were evaluated by TIDE, with a lower TIDE score suggesting lower possibility of immune escape and a more active immunotherapy response (Wang et al., 2022). Besides, the aneuploid score of TCGA samples was obtained from a published article (Huang & Fu, 2019), with higher aneuploidy score indicating higher possibility of evading immune surveillance and less active immunotherapy response.

GSEA

Pathway enrichment analysis was conducted using GSEA and the pathways with P_adj < 0.05 were defined as differential pathways. The HALLMARK pathway gene set was collected from the molecular feature database (https://www.gsea-msigdb.org/gsea/msigdb/) as a background gene set.

Identification of DEGs

The DEGs in S1 subtype and S2 subtype were selected with the limma package under |log2FoldChange(FC)| > log2(1.5) and P_adj < 0.05 (Song et al., 2023). The clusterProfiler R package (Lin et al., 2024) was employed in Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, and the pathways with P_adj < 0.05 were considered as significantly enriched (Song et al., 2023).

Development of a RiskScore model and validation

The number of prognostically relevant genes screened by performing univariate Cox proportional hazards regression (p < 0.05) were reduced by conducting LASSO Cox regression analysis in the R package “glmnet” (Li, Lu & Yin, 2022). Key genes and correlation coefficients were obtained through stepwise regression, and the formula of the RiskScore was as follow (Chu et al., 2023):

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum \mathrm{\beta }\mathrm{i}\ast \mathrm{E}\mathrm{x}\mathrm{P}\mathrm{i}.$$

βi and ExPi represent the coefficient of gene in Cox regression model and the gene expression value, respectively.

Based on the median RiskScore, low- and high-risk groups of COAD patients were classified. Then, the receiver operating characteristic (ROC) curve was developed with “timeROC” package (Lu et al., 2022). The RiskScore model was also validated in AC-ICAM cohort.

Construction of a nomogram and validation

Based on the results of univariate and multivariate Cox regression analyses, the prognostic nomogram was established by integrating Age, pathologic_M, pathologic_stage, and RiskScore using rms package (Liu et al., 2023). The accuracy and reliability of the nomogram were validated by the calibration curve and the decision curve analysis (DCA), respectively.

Genomic feature analysis

The data of SNV somatic mutation was processed using “maftools” R package (Mayakonda et al., 2018), and the differences in somatic gene mutation rates of the two different cohorts were compared using the mafCompare function.

Cell culture

Dulbecco’s modified Eagle medium (DEME) (Gibco-Thermo Fisher Scientific, Waltham, MA, USA) that supplemented with 10% fetal bovine serum (FBS), and 1% Penicillin and Streptomycin (BI, Hamagshimi, Israel) was used to culture human normal colonic epithelial cells (NCM460) obtained from Baidi Biotech Ltd. (Hangzhou, China) and human colon adenocarcinoma cells (SW1116) purchased from the Kalang biotechnology company (Shanghai, China) at 37 °C and 5% CO₂.

RT-qPCR

Cells of 80–90% growth were treated by trypsin and washed with PBS for collection. Then, the TRIzol reagent Kit (Vazyme, Nanjing, China) was added to isolate total RNA, which was reverse-transcribed and quantified using a commercial kit (Vazyme, Nanjing, China) for synthesizing cDNA based on the manufacturer’s specification. Finally, the qPCR reaction system was mixed by the SYBR Green PCR Master Mix (Invitrogen) for qRT-PCR on LightCycler 96 (Roche), and the ^ΔCT, ^ΔΔCT value, 2^−ΔΔCT value were calculated for each sample (Kang et al., 2022). The specific primers were shown in Table S2.

Wound healing and transwell assay

The si-CDKN2A (Forward: 5′-GCGCGGAGCCCAACTGCG-3′; Reverse: 5′-AGCGTCGTGCACGGGTCG-3′) obtained from the Sangon Biotech company (Shanghai, China) was transfected into cells by adding Lipofectamine 3000 (Invitrogen). The SW1116 cells with si-CDKN2A silencing was used for the wound healing assays, and 4 × 10⁶ cells were seeded into the 12-well plate containing DMEM. After reaching 95% confluence, the tip of a 100 μl pipette was employed to scratch cell layers, ensuring the same width of each wound. Next, the cells were rinsed in PBS reagent and transferred into complete medium with 1% FBS for further incubation at of 37 °C with 5% CO₂. The width of the cell injuries was measured with an inverted microscope immediately after the scratch and 48 hours (h) later.

For transwell test, the 12-well plates with 8-μm pore size chamber (Corning) were used. The cell suspension (200 μl) was supplemented into the upper chamber of the Matrigel-coated serum-free transwell, while the lower chamber was added with 500 μl of complete medium with 10% FBS. After incubating the transwell setup in a 5% CO₂ environment for 48 h at 37 °C, 0.1% crystal violet solution were added to dye the cells in the lower chamber for 15 minutes (min) and the number of cells were counted with an inverted microscope (Leica, Wetzlar, Germany).

Flow cytometry

The harvested cells were resuspended in PBS at a density of 1 × 10⁵ cells per 200 μL and subsequently stained with Annexin V-FITC and PI solution (C1062S, Beyotime, China) on ice for 30 min in the dark. Following washing the cells with PBS, BD FACSAria III flow cytometer (BD, Franklin Lakes, NJ, USA) was utilized to analyze the samples (Tian et al., 2023).

Statistical analysis

All statistical data were analyzed in R language (version 4.2.0 and Graph Prism 8.0). Differences between two-group continuous variables were calculated with Wilcoxon rank-sum test. The statistics of experimental data were performed by student’s t test. The correlation was calculated using Spearman method. Survival difference in different groups was compared by plotting Kaplan-Meier (KM) curves with the log-rank test and p < 0.05 represented statistically significant.

Two molecular subtypes of COAD were identified using PCD-related genes

Firstly, 21 PCD-related prognostic genes for COAD were obtained from TCGA and AC-ICAM cohorts (Figs. 1A and 1B). Based on consensus clustering analysis, 427 COAD samples in TCGA were divided into two subtypes (S1 and S2) when consensus matrix k = 2 (Fig. 1C). Specifically, S2 subtype showed a relatively favorable prognosis, while S1 subtype had a worse prognosis (Figs. 1D and 1E). Further analysis revealed that the protective genes (IL13, GLB1, IL13RA2, SLC39A8, FAS) were high-expressed in S2, while the risk genes were high-expressed in S1 (Figs. 1F and 1G). In addition, compared with S2 subtype, the distribution of clinical and pathological features in S1 and S2 subtypes in the TCGA dataset revealed that more patients in S1 subtype were in the pathologic_T3, pathologic_T4, pathologic_N1, pathologic_N2, pathologic_M1, pathologic_stage III, and pathologic_stage IV (Fig. 1H).

Differences in immune microenvironment and pathway characteristics between S1 and S2 subtypes

Next, the immune microenvironment between S1 subtype and S2 subtype was compared. CIBERSORT analysis showed that the abundance of T cell CD4+ memory, macrophage M1, and mast cell resting was high in S2, while S1 had more Tregs and macrophage M0 (Fig. 2A). MCP-counter analysis indicated that the immune cell scores of natural killer (NK) cell and cytotoxicity score in S2 were remarkably higher than those in S1 and the cancer-associated fibroblast and endothelial cell were high in S1 group (Fig. 2B). This indicated that the S2 subtype may have a stronger anti-tumor immune microenvironment, while the S1 subtype may have an immunosuppressive environment. GSEA showed that S1 subtype were predominantly enriched in pathways related to tumor invasion and progression, including myogenesis, angiogenesis, and epithelial mesenchymal transition (EMT) (Fig. 2C). In contrast, S2 subtype was notably enriched in the pathways linked with cell proliferation, including G2m Checkpoint, MYC Targets V1, and E2f Targets (Fig. 2D).

Identification of DEGs between S1 and S2 subtypes

DEGs analysis screened a total of 90 downregulated and 394 upregulated genes from S1 and S2 subtypes in the TCGA cohort (Fig. 3A). The heatmap displayed the top 50 DEGs (Fig. 3B). KEGG analysis showed that the upregulated genes were principally enriched in TGF−β signaling pathway, ECM−receptor interaction, and proteoglycans in cancer (Fig. 3C), while the downregulated expression genes were largely enriched in the chemokine signaling pathway, IL−17 signaling pathway, and cytokine−cytokine receptor interaction (Fig. 3D).

RiskScore model with a strong prognostic prediction performance

A total of 484 DEGs were subjected to univariate Cox proportional hazards regression. Then, the gene number was further narrowed using LASSO Cox and stepwise regression (Fig. 4A). Six key genes were obtained to establish a RiskScore model of “RiskScore=(0.208*HSPA1A) + (0.273*SEMA4C)+(−0.087*ATOH1)+(0.154*CDKN2A)+(0.163*ARHGAP4) + (−0.073*ZG16)”, with ATOH1 and ZG16 as the protective genes (hazard ratio < 1) and SEMA4C, HSPA1A, ARHGAP4, and CDKN2A as the risk genes (hazard ratio>1) (Figs. 4B and 4C). COAD patients were divided into high- and low-risk groups by the median RiskScore value. Furthermore, the ROC for the TCGA cohort showed that the RiskScore model had 1-, 2-, 3-, 4-, and 5-year AUC of 0.69, 0.7, 0.69, 0.72, and 0.69, respectively (Fig. 4D). In TCGA cohort, the low-risk group had better OS, disease-specific survival (DSS), and progression-free interval (PFI) in comparison to the high-risk group (Figs. 4E–4G). Moreover, the ROC analysis showed that the RiskScore model was highly robust in predicting 1- to 5-year survival in the AC-ICAM cohort, with AUC values of 0.72, 0.61, 0.6, 0.58, and 0.6 (Fig. 4H). Survival analysis in AC-ICAM cohort also showed a low OS rate and poor prognosis in high-risk group (Fig. 4I). Additionally, principal component analysis (PCA) revealed an obvious boundary between the high- and low-risk groups in the two cohorts (Fig. 4J). Differential expression analysis of the model genes demonstrated that ATOH1 and ZG16 were high-expressed protective factors in the low-risk group, while SEMA4C, HSPA1A, ARHGAP4, and CDKN2A were high-expressed risk factors in the high-risk group (Fig. 4J).

A nomogram was developed with a strong predictive performance

Analysis on the clinical grades in TCGA cohort showed that the RiskScore was higher in the pathologic_stage III and IV than the pathologic_stage I and II, higher in the pathologic_T3 and T4 than the pathologic_T1 and T2, higher in the pathologic_N1 and N2 than the pathologic_N0, higher in the pathologic_M1 than the pathologic_M0 (Fig. 5A). The distribution of clinical and pathological features between the two risk groups of COAD patients demonstrated that more patients were in the pathologic_T4, pathologic_N2, pathologic_M1, and pathologic_ stage IV in high-risk group (Fig. 5B). Moreover, univariate and multivariate Cox regression analysis on the clinical features and RiskScore confirmed that Age, the RiskScore, pathologic_M, pathologic_stage were independent factors for evaluating the prognosis of COAD (Figs. 5C and 5D). Then, a nomogram was developed. It was observed that RiskScore exhibited the greatest impact on the prediction of survival rate (Fig. 5E). Furthermore, the predicted calibration curves for 1-, 3-, and 5-year survival by the nomogram were close to the standard curves (Fig. 5F), indicating that the prediction by the nomogram was accurate. Additionally, the reliability of nomogram was assessed by DCA, which showed that the calibration curves indicative of the benefits of RiskScore and nomogram were notably higher than the extreme curves (Fig. 5G).

High-risk group exhibited strong immune escape ability and poor immunotherapy response

Compared to the high-risk group, CIBERSORT analysis showed that the abundance of mast cell resting, T cell CD4+ memory, and B cell Plasma was higher in low-risk group (Fig. 6A). MCP-counter analysis revealed the strongest positive correlation between the RiskScore and cancer-associated fibroblasts (CAFs) (Fig. 6B). Moreover, comparison on the expressions of seven immune checkpoint genes indicated that the expression of four immune checkpoint genes (PDCD1, HAVCR2, TIGIT, and LAG3) was higher in high-risk group (Fig. 6C), showing that high-risk COAD patients had strong immune escape ability and limited benefits from immune checkpoint inhibitors. Additionally, the RiskScore and TIDE score were positively correlated, with a correlation coefficient R of 0.39 (Fig. 6D). The number of respondents was higher in low-risk group (39.72%) than in high-risk group (26.29%) (Fig. 6E), showing a greater immunotherapy benefit to low-risk COAD patients. The aneuploid score was higher in high-risk group than that in low-risk group (Fig. 6F), which showed that tumor cells in high-risk patients may have greater genomic instability and may be more likely to evade immune surveillance and become resistant to treatment.

Differences in GSEA and somatic mutation between the two risk groups

Pathways such as EMT, angiogenesis, and interferon-γ response were activated in high-risk group but fatty acid metabolism was inhibited (Fig. 7A). Further analysis showed that TP53 exhibited a particularly higher somatic mutation rate in high-risk group than in low-risk group (Fig. 7B).

Expression of the model genes validated in vitro assay

The results of RT-qPCR showed that HSPA1A, SEMA4C, ARHGAP4 and CDKN2A were significantly overexpressed in the SW1116 cells, whereas ATOH1 and ZG16 were significantly downregulated in the cells (Fig. 8A). As CDKN2A was a high-risk gene, its expression level was significantly elevated in high-risk group, which showed a poor prognosis. Hence, CDKN2A was chosen as the subject of functional validation due to such a strong correlation. The wound healing and transwell assay revealed that silencing CDKN2A significantly impaired the tumor cell migration and invasion (Figs. 8B–8E). The apoptosis rate of cells was significantly elevated after CDKN2A knockdown in comparison to a low apoptosis rate in the control group, suggesting that CDKN2A may play a crucial role in the survival of colorectal cancer cells by inhibiting apoptosis (Fig. 8F). These results further supported the potential oncogenic function of CDKN2A in COAD progression.