Calculation of interaction-perturbation in the cell death interplay network
To generate the interaction-perturbation of all gene pairs from the cell death interplay network, we introduced a four-step pipeline as previously reported[10] (Fig. 1). CRC tissues from TCGA-CRC (n = 567) and normal colorectal tissues from GTEx database (n = 304) were regarded as tumor and normal input, respectively. The initial rank matrix was obtained based on the rank of each gene in a single sample, which was subsequently converted to the delta rank matrix using subtraction in the same direction of gene interactions. Furthermore, gene interaction was relatively stable and conservative in normal tissues[12, 18], thus serving as the benchmark for calculating the interaction-perturbation matrix in tumor tissues. Collectively, we generated a perturbation network consisting of 8403 edges and 1056 nodes. Previous evidence has demonstrated that biological networks are commonly scale-free distribution[38], which was in line with our constructed network (R =-0.963, P < 0.001; Figure S1A). Notably, relative to normal samples, tumor samples exhibited stronger perturbations (Figure S1B), and the range of interaction perturbation is also broader range, suggesting a higher variation (Figure S1C). These results indicated that the cell death interplay network was overall stable in normal samples, whereas it significantly perturbated in tumor samples.
Perturbation of cell death networks deciphered four heterogenous subtypes
Previous studies have demonstrated that network perturbation could more reliably better characterize the biological state of bulk tissues[10–12]. Here, representative perturbation features that can strikingly distinguish normal and tumoral tissues as well as maintain significant heterogeneity across tumor samples were retained for subtype discovery, which formed a sub-network including 4203 edges and 924 nodes. Specifically, this sub-network also fitted a scale-free distribution (R =-0.977, P < 0.001, Figure S1D). Afterward, Consensus clustering was performed with different k (k = 2–9) clusters according to the perturbation matrix of representative features. Results of cumulative distribution function (CDF) curve and consensus score matrix suggested that the optimal division was achieved when k = 4 (Fig. 2A-B). Thus, four subtypes were decoded from cell death network (CDN1, n = 128, 22.6%; CDN2, n = 170, 30%; CDN3, n = 170, 30%, and CDN4, n = 99, 17.4%). Kaplan-Meier analysis displayed significant survival differences among four subtypes, with CDN2 having the best prognosis and CDN4 possessing the worst prognosis with a 5-year survival rate of about 50% (P = 0.013, Fig. 2C).
To further validate the reproducibility and stability of four CDN subtypes in independent cross-platform cohorts, we retrieved three internal datasets to perform multiclass predictions using the NTP algorithm. Initially, 799 subtype-specific genes were identified via differential expression analysis (Table S5). Subsequently, cluster prediction was achieved by the NTP framework integrated 799 featured genes and three testing expression matrices. Encouragingly, the sample subtypes displayed similar transcriptome profiles across different cohorts (Figure S2A-C). In line with the prior findings, CDN4 presented dismal prognosis relative to other subtypes in all testing datasets (Fig. 2D-F). In addition to similar transcriptome and clinical traits, four subtypes also demonstrated analogical proportion (Figure S2D). Taken together, the above suggested that our CDN taxonomy was robust and reproductive in cross-platform cohorts.
Correlations of four subtypes with clinical characteristics and classical subtypes
Subsequently, we compare the distribution of clinical characteristics among four subtypes (Fig. 2G). CDN1 was closely associated with KRAS mutation (P < 0.001), whereas CDN2 was related to BRAF mutation (P = 0.004) and high level of microsatellite instability (MSI-H) (P < 0.001) (Figure S2F-G). In addition, CDN4 was related to advanced TNM stage (P < 0.05), consistent with the poor prognosis of CDN4.
To further compare our CDN taxonomy with previously developed CRC classifications, patients were reclassified according to criteria based on previous studies from Guinney et al.[5], Isella et al.[6], De Sousa E Melo et al.[7], Sadanandam et al.[8] and Marisa et al.[9]. Intriguingly, the results revealed strong associations between four CDN subtypes and previous classifications, indicating the molecular convergence in CRC (Fig. 2G). Specifically, CDN1 was linked to the canonical CMS2, CRCA4, CRCA4, CRIS-C/CRIS-D, and CIT1; CDN2 was associated with CMS1, CCS2, CRCA1, CRIS-A, and CIT2; CDN3 was enriched in CMS2, CCS1, CRCA3, CRIS-C, and CIT5. In contrast, CDN4 was related to more aggressive classification features, including CMS4, CCS3, CRCA5, CRIS-B, and CIT4. Of note, only a fraction (20%) of our signature genes overlapped with signature genes of previous CRC classifications, suggesting that our classification has specific biological characteristics as well as more room for exploration (Figure S2E).
Biological characteristics of four subtypes
To characterize the biological behaviors in four subtypes, we calculated the variation of 9997 pathways via the GSVA algorithm. The results indicated that biological functions differed dramatically among four subtypes (Fig. 3A). Specifically, CDN1 was endowed with moderate metabolic activity and low immune activity, which demonstrated significant enrichment in proliferation pathways, such as cell cycle, G2M checkpoint, and MYC targets. CDN2 was also enriched in proliferation pathways but displayed conspicuous enrichment in immune-activated and multiple cell death pathways. The metabolic pathways, such as lipid, vitamin, and galactose metabolism, were particularly evident in CDN3. CDN4 was associated with significant aggressive features, including cancer stem cell-related pathways, metastasis-related pathways, and multiple cancer-related signaling pathways. Subsequently, the enrichment of the four subtypes in 11 cell death pathways was explored (Fig. 3B), and we observed that CDN2 exhibited high levels in all cell death pathways, while CDN1 exhibited the lowest levels. CDN3 displayed moderate levels in all cell death-related pathways. CDN4 was mainly enriched in autophagy and lysosome-dependent cell death pathways. Overall, CDN1 was defined as a pronounced subtype because of the significant proliferative activity. CDN2 was defined as an immune subtype with strong cell death activity. CDN3 had high metabolic activity as metabolic activity. As a result of high cell stemness as well as high levels of invasive tumor pathways, CDN4 was defined as an aggressive subtype, along with high autophagic activity.
Multi-Comics landscape of four subtypes
To characterize the molecular profiles of four CDN subtypes at the mutational and CNV levels. First, we depicted the mutational landscape of four subtypes (Fig. 4A), which showed that CDN1 and CDN3 were prominent in APC, TP53, and KRAS mutations, while TTN, SYNE1, and MUC16 performed higher mutation frequencies in CDN2 and CDN4 (Fig. 4D). Notably, KRAS mutations can stimulate tumor cell proliferation and growth, which leads to a poor prognosis, while mutations are less abundant in CDN2. Overall, most driver genes in CDN2 have the highest mutation frequency, and tumor mutation burden (TMB) is also the highest of the four subtypes (Fig. 4A-B). Higher TMB tend to generate more tumor neoantigens (Fig. 4C), which are more likely to stimulate immune activation. In addition, the burden of CNV (gain or loss) at the level of bases, fragments, and chromosome arms was highest in CDN3 (Fig. 4E). That said, CDN3 was inclined to be CNV-driven, whereas CDN2 was inclined to be mutation-driven.
The immune landscape of four subtypes
Previous studies have demonstrated that the composition of immune cells within infiltrating CRC tumors is heterogeneous, and the key players are continuously subject to microenvironmental changes. Thus, to further explore the immune landscape of four subtypes, the ssGSEA was employed to measure the abundance infiltration of 28 immune cells. The results displayed that most immune cells infiltrated mainly in CDN2 and CDN4, with CDN3 as an intermediate location, while CDN1 had the lowest level of immune cell infiltration (Figure S3). Especially, activated effector cells that serve an instrumental role in antitumor immunity, such as activated CD8 T cells, activated CD4 T cells, and cytotoxic lymphocytes (Fig. 5A), were highly infiltrated in CDN2 and CDN4. In addition, immune checkpoints such as B7-CD28 family member proteins and TNF superfamily were highly expressed in CDN2 and CDN4 (Fig. 5B). Among them, the immune checkpoint CD40L binding to CD40 could increase the immunomodulatory ability of DCs, induce effector cell proliferation and enhance immunotoxicity to tumor cells, which further indicated that CDN2 and CDN4 had strong immune activation ability. Furthermore, the APS score also demonstrated the strong antigen-presenting ability of CDN2 and CDN4 (Fig. 5F). However, although CDN4 had a higher immune score, macrophages M2, regulatory T cells, and stromal cells were highly infiltrated in CDN4 with a higher stromal score, which demonstrated that CDN4 was an immune hot but suppressed tumor microenvironment. (Fig. 5C-E). CDN1 possesses highest tumor purity accompanied by low levels of immune infiltration (Fig. 5E). Moreover, CDN2 had higher levels of single nucleotide variant (SNV) neoantigen, indel neoantigen, and CDN4 was featured by higher TCR Richness, TCR Shannon, homologous recombination deficiency, intratumoral heterogeneity, and cancer-testis (CTA) scores (Fig. 5H). Taken together, CDN1 was defined as the immune desert subtype because of insufficient immune cell infiltration. CDN2 had high immune activity and was defined as the immune activating subtype. CDN3 is the intermediate immune subtype. With the strong immune activation but suppressed microenvironment, CDN4 was defined as the immune suppressive subtype.
Immunotherapy and potential drug targets
To systematically assess the immunotherapeutic potential of four CDN subtypes, we need to establish a comprehensive assessment of tumor immunity. Karasaki et al. proposed an immunogram for quantifying the cancer-immunity cycle[30]. Hence, we applied the immunogram to assess the potential for subtypes to benefit from immunotherapy (Fig. 5I). CDN1 and CDN3 exhibited low immune cycle scores, in line with previous findings indicating inadequate immune infiltration. CDN2 and CDN4 displayed higher levels of immune activation and antitumor related immune cycle scores, along with relatively high immune checkpoint expression. Thus, CDN2 and CDN4 are more likely to benefit from immunotherapy. In addition, we applied an investigational 18-gene signature termed TIS, which can be used to measure a pre-existing but suppressed adaptive immune response within tumors. As shown in Fig. 5G, CDN2 and CDN4 displayed higher levels of TIS scores. SubMap analysis also revealed significant expression similarity between CDN2 and responders in CAR-T, anti-PD-1, and anti-CTLA-4 treatment cohorts (Fig. 6A). Moreover, CDN2 and CDN4 might be sensitive to anti-PD-1/CTLA-4 combination immunotherapy (Fig. 6B). These results demonstrated that CDN2 might benefit from immunotherapy, while CDN4 with the poorest prognosis might benefit from multi-targeted immunotherapy. Overall, results indicated that CDN2 and CDN4 might derive clinical benefit from immunotherapy, whereas CDN1 and CDN3 may not fit it due to the high cost and underlying immune-related adverse events that accompany low response potential.
Previous results indicated that CDN1 and CDN3 do not benefit from immunotherapy. Therefore, additional therapeutic approaches are required to improve the clinical outcomes of CDN1 and CDN3. Based on the CTRP and PRISM databases, we employed the drug evaluation model to search for drugs with therapeutic potential in four subtypes. First, the robustness of drug response assessment needs to be validated. Four EGFR pathway inhibitors are known to have better clinical efficacy in patients with KRAS mutation[39], and the results are consistent with the clinical results through validation (Fig. 6C). Subsequently, 12 subtype-specific drugs were identified for CDN1, CDN2, and CDN4 (Fig. 6D-F), but CDN3 did not receive a corresponding subtype agent. Notably, we observed that some specific-subtype drugs inhibit tumors through cell death pathways. For example, dasatinib and YM-155 can both induce apoptosis to inhibit tumor growth. Furthermore, clofarabine, a potential drug of CDN2, can inhibit tumors by inducing autophagy and apoptosis, demonstrating the potential of multiple cell death pathways operating together in treating CRC.
AOC3: a predictive molecule of CDN4 indicate poor prognosis
To further explore the CDN4 subtype with the worst prognosis, a total of 30 CDN4 feature genes kept AUC > 0.9 in four independent prognostic cohorts (Table S6). The result of univariate Cox analysis based on CDN4 signature genes demonstrated that AOC3, CRYAB, PROCKLE1, and TSN1 were risk factors with prognostic significance and showed significant consistency in the four cohorts (Fig. 7A). However, only three genes were eligible for both conditions, and the mean C-index of AOC3 (mean C-index: 0.612, Fig. 7B) was the highest compared with CRYAB (mean C-index: 0.598) and TSN1 (mean C-index: 0.605). In cohorts with overall survival information from TCGA and GEO, the low AOC3 group performed a significantly better prognosis than the high AOC3 group (Fig. 7C-F). To further validate the prognostic potential of AOC3, we divided the samples into the high H-scores group and the low H-scores group based on the H-scores of AOC3 from 103 successfully stained cancer samples in human tissue microarrays (Fig. 7G). The survival analysis result was consistent with previous analyses that the high AOC3 group had a significantly worse prognosis than the low AOC3 group (Fig. 7H).