Ferroptosis subclasses in CRC
The flow chart of this study is shown in Figure 1A. A total of 461 colon adenocarcinoma (COAD) and 172 rectum adenocarcinomas (READ) patients from TCGA database and 114 CRC patients from the GEO database were finally enrolled. We merged TCGA and GEO data and adopted the ComBat method to remove the batch effect. Next, we extracted ferroptosis expression gene profiles and found that there were 59 ferroptosis-related genes expressed in CRC tissues. Distribution of FRGs by circosplot were shown in Figure 1B.
By applying consensus clustering on the gene expression profile of the above mentioned 59 FRGs on metadata set 747 CRC samples, four resulting clusters were defined as C1-C4(Figure2A, B). In addition, a high similarity of gene expression pattern was observed within each subgroup based on consensus matrix k=4(Figure2C). Survival analysis indicated a significant prognostic difference between four subgroups in metadate set (Figure 2D, p=0.034). Next, we developed a heatmap demonstrating the relationship between stage and 59 genes (Figure 2E). To validate the subclasses’ assignments, we performed the correlation analysis of 59 genes, the results showed that there were co-expression correlations among some FRGs (Figure 3A). Principal component analysis (PCA) was applied to decrease the dimension of features and we found the subtype designations were largely concordant with two-dimensional PCA distribution patterns (Figure3B). Hence, it was reasonable to classify FRGs into four subclasses.
Transcriptomes of the CRC subclasses
To further characterize the differentially expressed FRGs (DEGs)in different subtypes, ANOVA differential expression analysis were performed. Gene expression differences were considered significant if the adjusted P value was < 0.05('BH' correction), the results indicated that all 59 FRGs were differentially expressed in all of these subtypes (Figure4).
To elucidate the biological functions and pathways of the 59 DEGs in CRC. Gene Oncology enrichment analysis of the signature genes was conducted using the CLUSTERPROFILER package, and significantly enriched biological processes are shown in Figure 5. GO enrichment analysis showed that these FRGs were significantly associated with biological processes such as metabolic processes and oxidative stress, and KEGG pathway analysis revealed that they were associated with metabolic pathways such as iron-dependent cell death and glutathione metabolism.
Tumor immune microenvironment characteristics among four subclasses
To uncover the immune heterogeneity among the four subtypes, we resorted to immune-related tools. Through the CIBERSORT method, we evaluated the differences in the immune infiltration of 22 immune cell types between four subclasses in CRC patients. As illustrated in Figure 6, the infiltration of 18 immune cells were different between the four subclasses. Compared with the other three subclasses, tumor in C3 had higher abundance of 7 immune cell populations (B cells, M2 macrophages, resting mast cells, monocytes, NK cells, plasma cells, CD8+T cells). In addition, C3 also exhibited low infiltration of M1 macrophages, activated mast cells, neutrophils, Tfh cells and Treg cells (Figure 6). We furtherly investigated stromal and immune scores based on the ESTIMATE method. The results indicated that tumor in C3 had higher immune and stromal scores but lower tumor purity (Figure7).
Assessment of immunotherapeutic response in subclasses
Immunotherapy have emerged as a promising strategy for the treatment of many cancers. We subsequently investigated the expression of PD-L1 which is a pivotal immunomodulator. As shown in Figure 8A, the expression of PD-L1 were different between the four subclasses, and tumor in C3 had higher expression of PD-L1 expression than other subclasses. In addition, we investigated the correlation between PD-L1 expression and FRGs in CRC patients. The results revealed that the expression of PD-L1 was significantly correlated with multiple FRGs such as ACSL4, ALOX5, HOMX1 and IREB2(Figure8B, FigureS1). Cancer stem cells are known to play a critical role in the growth, metastasis, and recurrence of tumors. Recent studies have reported that higher stemness indices could limit the antitumor immune responses and be associated with decreased PD-L1 expression[28]. We herein applied an innovative one-class logistic regression (OCLR) machine-learning algorithm to calculate stemness indices of subclasses. As shown in Figure9, there was a significant difference in the stemness index among the subtypes, and we observed that patients in C3 had lower stemness indices value than other subtypes. Thus, these data support our analysis that patients in C3 may have better prognosis, and might be more promising in response to immunotherapy.
Difference sensitivity to antitumor drugs for CRC subclasses
By applying the R package pRRophetic, we further explored the association between the CRC subclasses and sensitivity toward antitumor drugs. For gemcitabine, there were significant differences between the four subtypes(p<0.0001). For eriotinib, C3 exhibited significant association of eriotinib-resistant group. For bortezomib, the results showed significant differences between some subtypes (Figure10). In addition, other drugs such as cytarabine, dasatinib and sorafenib were explored, and significant differences were observed between some subtypes (Supplementary Figure). These data demonstrated the CRC subclasses might provide precise therapeutic approaches in the further.
Correlation of the CRC subclasses with mutations and CNV differences
The tumoral genomic landscape has been proven to be associated with the application and efficacy of targeted therapies and immunotherapy. We further evaluated the difference of gene mutations and copy number variation (CNV) among the four subtypes. The results revealed that distinct subtypes had different subclasses tend to have different mutational characteristics. For example,49% samples in C1 had TP53 mutation while 62% in C2 and 71% in C4(Figure11A). What’s more, we further obtained gistic scores of different subtypes according to the genepattern website, and the results indicated that there was a significant difference in mutation frequency between the four subtypes. These data could protect the samples in different clusters from the option of individualized therapeutic agents.
Correlation of CRC subclasses with clinical parameters
In order to elucidate whether the ferroptosis-related subclasses for CRC is independent indicator when considering other conventional clinical parameters. Based on the integrated clinical information downloaded from COAD and READ in TCGA with the clinical information obtained from GEO, we obtained the clinical information of the metadate dataset. By applying the RCIRCOS method, we compared the distribution of ferroptosis subclasses with the patients’ stage, age, gender, tissue or organ of origin, and BMI. As shown in Figure 10, ferroptosis subclasses do not show significant crossover with other clinical parameters, and have significantly difference with the parameters. Hence, the subclass may serve as a potential molecular subtype which can provide an important prognostic and therapeutic value for CRC.