Graphical abstract graphic file with name fx1.jpg [57]Open in a new tab Highlights * • Proteomic subtypes reveal immune regulation associated with drug sensitivity * • HSF1 can increase DNA damage repair, thus inducing resistance to radiation * • HDAC6 regulated by NFYC can synergistically induce cetuximab resistance * • Prognostic models with high accuracy are constructed for response prediction __________________________________________________________________ Molecular characterization associated with therapeutic response in CRC are incompletely characterized. Li et al. identify proteomic subtypes associated with distinct biological features and therapeutic responses, explore the potential mechanisms related to sensitivity and resistance, and establish models for response prediction to improve effectiveness of clinical treatment. Introduction Colorectal cancer (CRC) is one of the most common malignant tumors of the digestive system and a leading cause of cancer-related death worldwide. Unlike other cancers, no single risk factor accounts for most cases of CRC.[58]^1 Preoperative chemotherapy (CT; mainly FOLFOX [folinic acid, fluorouracil, and oxaliplatin] therapy), chemoradiation therapy (mainly FOLFOX plus radiation), and CT plus targeted therapy (mainly FOLFOX plus cetuximab or bevacizumab) followed by surgery are currently the recommended standard treatments for locally advanced CRC. However, only half of patients with CRC respond to FOLFOX CT,[59]^2 and approximately 10%–30% have pathologic complete response after receiving chemoradiation therapy.[60]^3 In recent years, drug development has focused on agents targeting proteins expressed strongly or exclusively by tumor cells. Preclinical and clinical studies have shown that cetuximab and bevacizumab are the approved antibodies separately targeting EGFR (epidermal growth factor receptor) and VEGF (vascular endothelial growth factor). According to the pan-Asian adapted European Society for Medical Oncology (ESMO) consensus recommendations, the testing of RAS mutational status should be carried out on the tumors of patients with CRC at diagnosis.[61]^4 As previous trials demonstrated, cetuximab plus CT has been typically used in the first-line treatment of RAS wild-type metastatic CRC.[62]^5^,[63]^6^,[64]^7^,[65]^8 A randomized controlled trial revealed that cetuximab plus CT improved the objective response rate (ORR) (57.1% versus 29.4%) compared with CT alone for patients with initially unresectable KRAS wild-type colorectal liver metastases.[66]^8 The BECOME randomized controlled trial showed that bevacizumab plus mFOLFOX6 improved response rates (54.5% versus 36.7%) compared with mFOLFOX6 alone for patients with initially unresectable RAS-mutant colorectal liver metastases.[67]^9 Drug resistance occurs extensively in different therapies. Thus, stratification of patients with CRC for appropriate therapy is indispensable for better therapeutic results. Several molecular classifications have been proposed to connect molecular patterns to clinical features. The Cancer Genome Atlas (TCGA) project mapped the genomic landscape of CRC and revealed that three-quarters of hypermutated CRC patients had high microsatellite instability combined with hypermethylation and MLH1 silencing.[68]^10 As proteins link genotypes to phenotypes, the Clinical Proteomic Tumor Analysis Consortium (CPTAC) presented the first integrated proteogenomic characterization of the TCGA CRC cohort and identified five proteomic subtypes featured with distinct mutation, methylation, and protein expression patterns associated with different clinical outcomes.[69]^11 Next, a collaborative paradigm that combined six independent classification systems was identified as four consensus molecular subtypes with distinguishing features.[70]^12 Recently, proteogenomic integration studies also suggested potential features related to drug responses.[71]^13^,[72]^14 Overall, these proteogenomic studies could provide approaches to the critical need for biomarkers to predict therapy response among CRC patients. Herein, we constructed a comprehensive CRC cohort, including a discovery cohort (N = 124) and a validation cohort (N = 130), which covered three clinical therapy subcohorts, including a CT subcohort, a chemoradiotherapy (CRT) subcohort, and a CT combined with targeted therapy (CTT) subcohort for a systematic identification of potential therapeutic opportunities. Proteomic clustering resulted in four CRC molecular subtypes with distinct biological functions and therapeutic responses, which showed high concordance with TCGA proteomic subtypes. We further performed comparative analysis between sensitive (S) and non-sensitive (NS) groups separately in different therapy subcohorts and explored the different resistance/sensitivity mechanisms. Finally, we developed predictive models with high accuracy to predict therapeutic response, which were validated using parallel reaction monitoring (PRM) assay in an independent validation cohort composed of 130 CRC patients receiving FOLFOX therapy (N = 45), FOLFOX plus radiotherapy (N = 40), or FOLFOX plus cetuximab (N = 45). Collectively, this study provides comprehensive molecular subtyping associated with therapeutic response related to distinct therapies and implicates the underlying regulatory mechanisms of CT, CRT, and targeted therapies that may benefit clinical practice. Results The CRC cohort: The response to chemoradiation and targeted therapies in patients with CRC To investigate the proteomic patterns associated with the response to CT and targeted therapies, we assembled a cohort composed of 254 CRC patients, including a discovery cohort (N = 124) and a validation cohort (N = 130), who received CT (mainly FOLFOX therapy), CRT (mainly FOLFOX plus radiotherapy), or CTT (mainly FOLFOX plus cetuximab) ([73]Figure S1A). The detailed therapy regimens and response evaluation are shown in the [74]STAR Methods. In the discovery cohort, the CT subcohort (N = 48) was grouped into a CT-S group (CT[S]; N = 19) and a CT-NS group (CT[NS]; N = 29); the CRT subcohort (N = 45) was grouped into a CRT-S group (CRT[S]; N = 27) and a CRT-NS group (CRT[NS]; N = 18); and the CTT subcohort (N = 31) was grouped into a CTT-S group (CTT[S]; N = 18) and a CTT-NS group (CTT[NS]; N = 13) ([75]Figure S1B). For the independent validation cohort, 130 CRC treatment-naive patients’ formalin-fixed paraffin-embedded (FFPE) samples were collected, including FOLFOX (N = 45; N [S] = 18, N [NS] = 27), FOLFOX plus radiotherapy (N = 40; N [S] = 21, N [NS] = 19), and FOLFOX plus cetuximab (N = 45; N [S] = 23, N [NS] = 22) therapies. Archival FFPE tissues were collected from 254 therapy-naive patients with CRC, then subsequently analyzed using a mass spectrometry (MS)-based, label-free, data-dependent acquisition (DDA) quantitative proteomics method in the discovery cohort[76]^15^,[77]^16 and a PRM targeted quantitative proteomics method in the independent validation cohort ([78]Figure 1A).[79]^17^,[80]^18 Formal assessment showed no baseline differences in gender and age among the CT, CRT, and CTT subcohorts. With respect to primary site, we found in the CRT subcohort that 91.1% patients had rectal cancer and were receiving CRT, and the remaining 8.9% had colon cancer ([81]Figure 1B; [82]Table S1). Figure 1. [83]Figure 1 [84]Open in a new tab Summary of the proteomic analysis of colorectal cancer (A) The proteomics workflow involved three modules: cohort construction (including the discovery cohort and validation cohort), proteomic profiling, and data analysis. (B) The clinical characteristic of the three subcohorts. (C) Protein identification and protein overlap of the three subcohorts. Data are shown as mean ± SD. (D) PCA of the three subcohorts. (E) Dynamic range of the protein identification of each sample according to the descending sort of protein abundance. The highest and lowest abundance proteins and some known biomarkers associated with CRC therapy response are shown. (F) Protein identification and protein overlap of S and NS groups. (G) Heatmap showing no difference of these proteins’ expression between the S and NS groups. Proteomic analysis of the CRC cohort The MS platform was stable and repeatable as judged by quality control runs during the entire data collection period ([85]Figure S1C; [86]STAR Methods). Proteomics measurement resulted in a total of 11,558 gene products (GPs) in all samples and 6,519–7,372 GPs in each sample ([87]Figures 1C–1F, [88]S1F, and S1G; [89]Table S1). Among three subcohorts, there was no difference in the proteomic coverage ([90]Figure 1C), and principal-component analysis (PCA) also revealed no batch effects ([91]Figures 1D and [92]S1H). Some known biomarkers were reported to be involved mainly in tumor progression and had implications for the clinical therapy strategies of CRC, such as EZR,[93]^19 DPYD,[94]^20^,[95]^21 TYMP,[96]^14^,[97]^22 KRAS, NRAS,[98]^23 PTEN, BRAF, PIK3CA,[99]^23^,[100]^24 and TP53.[101]^25^,[102]^26 In our study, we found that these reported biomarkers were highly identified in this CRC cohort ([103]Figure 1E). Furthermore, differential analysis demonstrated that TP53 and TYMP showed higher increase in the tumor samples compared with the normal samples both in the CPTAC and TCGA cohorts, validating the significance of the CRC biomarkers and indicating their potential association with therapy response ([104]Figure S1I). Then, we preliminarily explored whether these biomarkers were exactly associated with therapy response among CRC patients on the basis of our proteome data. After ensuring no major differences in the proteomic coverage between the S (10,599 GPs) and NS (10,616 GPs) groups ([105]Figure 1F), we further compared these biomarkers between the S and NS groups and found no significant difference in our cohort ([106]Figure 1G), indicating the unmet need to search for potential biomarkers for CRC therapy. Our study has so far presented a comprehensive view of the proteomic landscape of this CRC cohort associated with first-line therapies, aiming to find potential biomarkers for distinguishing S patients from NS patients. Proteomic subtyping of the CRC cohort and its association with therapy response Proteins have been regarded as the “executors of life,” but we have no direct evidence that the proteome pattern can distinguish different therapy responses. Herein, we first applied consensus clustering analysis on the discovery cohort and identified four proteomic subtypes, in which 57, 48, 9, and 9 patients were grouped into subtypes G-I, G-II, G-III, and G-IV, respectively ([107]Figures 2A and [108]S2A; [109]Table S2). The proteomic subtypes displayed distinct clinical characteristics, especially in therapy response (p < 0.05, Fisher’s exact test) ([110]Figures 2B and [111]S2B). The G-III subtype exhibited a higher proportion (77.8%) of S patients, while the G-IV subtype exhibited a higher proportion (88.9%) of NS patients ([112]Figure 2C). Comparison of this subtype with TCGA subtype revealed high concordance, and G-IV subtype conferred the worst prognosis ([113]Figures S2C and S2D). Further functional enrichment analysis determined the dominant pathways of each subtype ([114]Figure 2D). The G-III subtype, prone to sensitivity to therapy, was featured mainly by immune-related pathways, while the G-IV subtype, prone to non-sensitivity to therapy, was dominant in the RAS-MEK-ERK signaling pathway ([115]Figures 2D and 2E). The representative differentially expressed proteins (DEPs) with implications in these overrepresented pathways are shown in [116]Figure 2E. These findings provided the direct basis for the association of proteome pattern and therapy response, which prompted us to find more relevant molecular features for response. Figure 2. [117]Figure 2 [118]Open in a new tab Proteomic subtyping of the colorectal cancer cohort and its association with clinical characteristics (A) The relative abundance of the signature proteins and biological functions of four subtypes. (B) The association of four proteomic subtypes with clinical characteristics. (C) Distributions of S and NS groups among four proteomic subtypes. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of four subtypes. (E) Differentially expressed signatures and signaling cascades involved in G-III and G-IV. (F) GSEA of complement cascade and neutrophil degranulation between the S and NS groups. (G) Signature proteins up-regulated in G-III subtype and in the S group and its validation of the association with drug response in DepMap database. (H) Differential expression of these proteins on the basis of PRM quantification in an independent cohort. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001. To validate these biological pathways associated with therapy response, we performed a direct comparative analysis between the S and NS groups in the discovery cohort. As a result, the top enriched gene set enrichment analysis (GSEA) pathways in the S group included complement cascade and neutrophil degranulation ([119]Figure 2F). Consistently, these results revealed the association of immune-related pathways with therapy sensitivity. Importantly, among the 1,663 proteins up-regulated in G-III subtype, 573 proteins showed significant increase in the S group. To investigate the association of proteins expression with therapy sensitivity, we explored the drug susceptibilities of these 573 proteins using CRC cell line data from the Cancer Dependency Map Project (DepMap). The results demonstrated that 219 proteins showed associations with drug sensitivity of 5-fluorouracil (5-FU) and oxaliplatin, respectively. Of them, 106 overlapping proteins were validated and associated with drug sensitivity of both 5-FU and oxaliplatin. Furthermore, we found a group of proteins involved in neutrophil degradation (RNASE2, PTPRJ, and CEACAM8) and complement/coagulation cascades (SERPINA5, C7, and F10), and the abundance of these proteins showed negative correlation with sensitivity of drugs (5-FU and oxaliplatin) (Pearson’s r < 0), indicating that patients featured with immune response pathway were more likely to benefit from CT or targeted therapies ([120]Figure 2G). These results demonstrated that the immune bioprocess contributed to therapy response of CRC patients. To further validate the association between immune-regulated proteins and sensitivity of drugs, we used the PRM assay, a powerful targeted MS approach adopted in recent proteomic research,[121]^17^,[122]^18 to quantify these proteins in the independent validation cohort (including 62 S patients and 68 NS patients). To validate the technically reproductivity, we also performed pairwise proteomic measurement of PRM and DDA. Interestingly, we observed a relative high correlation, with an average Pearson correlation of 0.90 between PRM and DDA ([123]Figure S2E). In the PRM analysis ([124]STAR Methods), we observed that these proteins were more highly expressed, with fold change (FC) > 2 in the S group compared with the NS group (p < 0.05, Wilcoxon rank-sum test) ([125]Figure 2H; [126]Table S2), validating the association of immune-related proteins with drug sensitivity. Overall, the proteomic subtypes provided a thorough understanding of different therapies, which was crucial to develop an optimized multimodality treatment plan. PLA2G4A-mediated cross-talk signaling associated with chemoresistance of CRC featured with RAS mutation To further explore the mechanism of sensitivity/resistance of specific treatments, we performed a comparative proteomic analysis on the therapy-naive samples from the S (N = 19) and NS (N = 29) patients included in the CT subcohort (mainly FOLFOX CT). The main clinical characteristics and RAS mutation in CT subcohort are shown in [127]Figure 3A. We found 169 DEPs in the CT subcohort (p < 0.05, Wilcoxon test; Benjamini-Hochberg [BH]-adjusted p < 0.05) ([128]Figure S3A). For patients receiving FOLFOX CT included in the CT subcohort (79% [N = 38]), we identified 137 and 58 DEPs in the S (N = 16) and NS (N = 22) groups, respectively ([129]Figure S3B). Because 79% of patients received FOLFOX therapy, accounting for the greatest proportion in the CT subcohort, ConsensusPathDB pathway enrichment analysis revealed consistent results that the T cell receptor (TCR) signaling, innate immune system, FcγR-mediated phagocytosis, chemokine signaling, and autophagy pathways were significantly enriched in the S group, while in the NS group, the overrepresented pathways were glycerophospholipid metabolism and purine/pyrimidine metabolism ([130]Figure 3B; [131]Table S3). In summary, immune regulation signaling pathways were associated with CT sensitivity, while metabolism bioprocesses were associated with CT resistance. Figure 3. [132]Figure 3 [133]Open in a new tab Differential expression of proteins and signaling pathways in the sensitive and non-sensitive groups of the chemotherapy subcohort (A) Detailed clinical characteristics of chemotherapy subcohort. (B) Pathway enrichment results of chemotherapy subcohort and FOLFOX subcohort. (C) Differential expression of proteins involved in pyrimidine/purine metabolism and glycerophospholipid metabolism. (D) Differential expression of the S/NS-sig quantified using the PRM approach (S, N = 18; NS, N = 27). (E) Volcano plot showing the correlation of these proteins and sensitivity of drugs. (F) Distribution of S/NS groups between RAS mutation group and RAS wild-type group. (G) Up-regulated pathways and related proteins in the RAS mutation group. (H) Differential expression of PLA2G4A in the RAS mutation/wild-type groups in our cohort (left) (RAS mutation, N = 39; RAS wild-type, N = 53) and in CRC cell lines (right) (n = 4 biologically independent samples in each group, respectively). (I) Differential expression of PLA2G4A in the S (N = 16) and NS (N = 23) groups of FOLFOX for CRC patients with RAS mutation. p values were calculated using two-sided Wilcoxon rank-sum tests (C, D, H, and I); ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001. Boxplots show median (central line), upper and lower quartiles (box limits), and 1.5 × interquartile range (whiskers) (D, H, and I). Purine and pyrimidine metabolism plays a critical role in the proliferation of cancer cells[134]^27; thus, impaired purine metabolism is associated with cancer progression. The key enzymes involved in purine/pyrimidine metabolism were up-regulated in the NS group in CRC patients receiving FOLFOX CT ([135]Figure 3C). Consistent with a previous study,[136]^28 resistance to 5-FU was linked to elevated expression of the main target TYMS, which catalyzed the de novo pathway for the production of deoxythymidine monophosphate. In addition, we found that proteins involved in glycerophospholipid metabolism were significantly up-regulated in the NS group. PRM analysis further demonstrated that the proteins involved in the purine/pyrimidine metabolism (POLE, POLR1D, DCK, and TYMS) and glycerophospholipid metabolism (GP1D, PCYT2, and PLA2G4A) pathways were also significantly up-regulated in the NS group in the independent cohort composed of 45 CRC patients receiving FOLFOX CT (FC [NS/S] > 2; p < 0.05, Wilcoxon rank-sum test) ([137]Figure 3D; [138]Table S3). Moreover, we evaluated the drug susceptibility of these proteins in response to 5-FU and oxaliplatin, which were the main components of FOLFOX regimens, in CRC cell line data from DepMap. The correlation analysis showed that PLA2G4A, a member of the cytosolic phospholipase A2 group IV family, showed the most significant positive correlation with sensitivity of both 5-FU and oxaliplatin (p < 0.05, r > 0) ([139]Figure 3E). As reported, the expression of PLA2G4A was possibly mediated by RAS mutation.[140]^29 The three RAS GTPases (KRAS, NRAS, and HRAS), particularly the family member KRAS, are the most commonly mutated oncogenes in cancer.[141]^30 As reported, the distribution frequencies of gene mutations in CRCs were ∼40% KRAS and 3%–5% NRAS, respectively.[142]^31 In this study, we found that 73% CRC patients with RAS mutations showed resistance to FOLFOX CT, indicating that patients with RAS mutations were unlikely benefit from FOLFOX CT ([143]Figure 3F). To investigate the downstream pathways regulated by RAS mutation, we then performed comparative proteomic analysis between patients with RAS mutation and wild-type. We found 289 proteins significantly elevated in the RAS mutation group, which were dominantly enriched in the MAPK signaling, FcγR-mediated phagocytosis, VEGF signaling, Ras signaling, PI3K-Akt signaling, and glycerophospholipid metabolism pathways (p < 0.05). Among these up-regulated proteins and overrepresented pathways, we found that PLA2G4A was involved in the most pathway items ([144]Figure 3G). We observed higher expression of PLA2G4A in CRC patients with RAS mutation (p < 0.05), which was also validated in other CRC cohort ([145]Figure 3H). Furthermore, differential analysis revealed RAS mutated NS patients had a 6.8-fold increase of PLA2G4A expression ([146]Figure 3I). Overall, we propose that besides purine and pyrimidine metabolism, PLA2G4A-mediated glycerophospholipid metabolism had a significant association with chemoresistance of CRC patients, especially for CRC patients with RAS mutations. These findings provided us the potential therapeutic strategies by targeting purine/pyrimidine and glycerophospholipid metabolism, which might be used to reprogram cancer metabolism. HSF1 regulated one-carbon metabolism and cell cycle and was associated with radiotherapy resistance of rectal cancer Besides CT, CRT is also advocated as the standard treatment for locally advanced rectal cancer. In this study, the CRT subcohort included 45 patients (N [S] = 27, N [NS] = 18), 66.7% of whom received FOLFOX plus radiotherapy (N [S] = 19, N [NS] = 11). We performed a comparative proteomic analysis on S and NS patients and identified 266 DEPs and 310 DEPs, respectively, in the CRT subcohort and the FOLFOX plus radiotherapy subcohort (p < 0.05, BH-adjusted p < 0.05) ([147]Figures S3C and S3D; [148]Table S4). The ConsensusPathDB pathway enrichment analysis revealed immune-related pathways, such as neutrophil degranulation, chemokine signaling, focal adhesion, ECM, and platelet activation, were significantly enriched in the S group (p < 0.05); while for the NS group, the up-regulated proteins were mainly involved in the metabolism-related pathways, such as amino acid metabolism, membrane trafficking, vesicle-mediated transport, and metabolism of vitamins and cofactors. We found consistent pathway enrichment results in the FOLFOX plus radiotherapy subcohort ([149]Figures 4A and 4B; [150]Table S4), and the differential expression of the proteins enriched in these pathways is shown in [151]Figure 4C. Figure 4. [152]Figure 4 [153]Open in a new tab Differential expression of proteins and signaling pathways in the sensitive and non-sensitive groups of the chemoradiotherapy subcohort (A) The pathway enrichment results are shown. (B) The proportion of patients treated with FOLFOX plus radiotherapy in the chemoradiotherapy subcohort and the consistent pathway enrichment in S and NS groups. (C) Differential expression of proteins involved in the consistent featured pathways. (D) The interaction network of proteins involved in the one-carbon metabolism pathway. (E) Differential expression of TFs in the FOLFOX plus radiotherapy subcohort. (F) Differential expression of HSF1 in the S and NS groups in the [154]GSE68204 dataset (S, N = 18; NS, N = 17). Boxplots show median (central line), upper and lower quartiles (box limits), and 1.5 × interquartile range (whiskers). (G) Spearman correlation of HSF1 and one-carbon metabolism ssGSEA score. (H) ROC curve of the combined predictive model signatures. (I) Differential expression of the S/NS-sig quantified using the PRM approach in an independent cohort (two-sided Wilcoxon rank-sum test; S, N = 21; NS, N = 19). (J) Qualification of HSF1 stained by IHC in representative examples. The scale bar indicates 100 μm. Data are analyzed using two-sided Student’s t test and shown as mean ± SD (n = 3 independent experiments). (K) Schematic illustration indicating the potential regulation HSF1 in response to chemoradiotherapy. Cancer metabolism is usually characterized by altered metabolic demand, nutritional supply, and metabolic enzyme regulation; targeting the specificity of cancer metabolism becomes a major strategy in cancer treatment. The proteins involved in metabolism of vitamins and cofactors play the important effect on the regulation of folate-mediated one-carbon metabolism, providing the necessary nucleotides and one carbon group for DNA replication. We found that proteins involved in one-carbon metabolism, including AHCY, ALDH1L1, ALDH1L2, GLRX, MTHFD1L, MTHFD2, and SHMT2, were significantly increased in NS patients receiving FOLFOX plus radiotherapy ([155]Figure 4D). Then we investigated the transcription regulation difference associated with CRT response and found that 8 transcription factors (TFs) were significantly different between the S and NS groups (FC > 2; p < 0.05, Wilcoxon rank-sum test) ([156]Figure 4E). Among these TFs, HSF1 was identified with the most significant difference. In addition, the association of HSF1 with radiotherapy was further validated in an independent published cohort composed of 59 rectal cancer patients receiving radiotherapy (32 NS patients and 27 S patients) on the basis of transcriptional data from the [157]GSE68204 dataset ([158]Figure 4F) ([159]STAR Methods).[160]^32 Moreover, HSF1 showed the highest correlation with one-carbon metabolism single-sample GSEA (ssGSEA) score (r = 0.42, p = 0.02) ([161]Figure 4G). A combined predictive model of these signatures involved in one-carbon metabolism regulation, including AHCY, ALDH1L1, ALDH1L2, MTHFD1L, SHMT2, GLRX, MTHFD2, and HSF1, showed good performance (area under the curve [AUC] = 0.89) in distinguishing S patients from NS patients receiving FOLFOX plus radiotherapy ([162]Figure 4H). In the CRT subcohort, there were 3 patients with colon cancer receiving CRT, which was the atypical case. To be more rigorous, we performed a much more relevant analysis on the patients with rectal cancer. Consistently, we found a similar pathway enrichment in the S and NS groups, in the rectal cancer patients included in the CRT subcohort and the FOLFOX plus radiotherapy subcohort (after removing the 3 non-rectal-cancer patients) ([163]Figure S3E). The difference of HSF1 and performance of the predictive model were also validated ([164]Figures S3F and S3G). Overall, our study demonstrated good consistency of these findings related to CRT in the rectal cancer patients. More important, these key proteins were further validated using the PRM approach in the independent cohort composed of 40 rectal cancer patients receiving FOLFOX plus radiotherapy (FC [NS/S] > 2; p < 0.05, Wilcoxon rank-sum test) ([165]Figure 4I). Immunohistochemistry (IHC) staining showed a significant difference of HSF1 between the S and NS groups in patients receiving CRT (p < 0.0001, Student’s t test) ([166]Figure 4J), which verified the aforementioned observations that HSF1 expression negatively responded to radiotherapy. Overall, these results indicated that HSF1 could transcriptionally regulate one-carbon metabolism and was associated with resistance to CRT ([167]Figure 4K). To explore the potential regulation of HSF1 in response to radiotherapy, we performed further in vitro validation experiment to confirm the association of HSF1 with CRT resistance. We first constructed an HSF1-overexpressing stable cell line (HSF1-OE) in SW480 and SW620 and an HSF1-knockdown cell line (HSF1-KD) in HCT116 and RKO ([168]Figures S4A and S4B). Compared with the control group, overexpression of HSF1 promoted the proliferation of SW480 and SW620 ([169]Figure 5A), while we observed an inhibition of survival in the HSF1-KD group ([170]Figure S4C), which was consistent with previous findings.[171]^33 We then determined the phenotype of radioresistance of CRC cells with differential expression of HSF1. Cell viability assay showed that the HSF1-OE group was more resistant to radiation compared with the control group, resulting in a significant increase in survival ([172]Figure 5B). As reported, phospho-H2AX (γH2AX) levels and foci are sensitive markers for DNA damage.[173]^34 Therefore, we further investigated the influence of HSF1 on γH2AX levels and foci in X-ray-treated CRC cells (SW480 and SW620) by immunofluorescence assay. The results showed there was no difference of foci formation after 2 h between the HSF1-OE and control groups with X-ray treatment, which demonstrated the same initial dose of radiation, while after 24 h, we observed lower γH2AX levels and fewer foci in the HSF1-OE group than in the control group at 24 h posttreatment with 2 Gy X-rays ([174]Figure 5C). Overall, HSF1-overexpressing cells showed more radioresistance, accompanied by increased DNA repair after radiation. Figure 5. [175]Figure 5 [176]Open in a new tab Effect of HSF1 expression on response to radiotherapy in CRC cell lines (A) Effect of HSF1 overexpression on proliferation in CRC cells (n = 3 independent experiments, two-sided Student’s t test, mean ± SD). (B) CRC cells were exposed to increasing doses of ionizing radiation (0, 2, 4, and 8 Gy), and surviving cells were quantified by CCK8 (n = 6 independent experiments, two-sided Student’s t test, mean ± SD). (C) Representative γH2AX foci present in HSF1-OE and control cells. Differential proteins and KEGG pathway alteration of the pre-irradiation (D) and post-irradiation (E) HSF1-OE group and HSF1-KD group. (F) Differential proteins and KEGG pathway alteration of post-irradiation 1, 2, and 12 h compared with pre-irradiation HSF1-OE group. (G and H) Boxplot showing the cell cycle and autophagy ssGSEA scores of the HSF1-OE/KD groups pre-irradiation and 1, 2, and 12 h post-irradiation (ANOVA, n = 3 biologically independent samples in each cell line). Boxplots show median (central line), upper and lower quartiles (box limits), and 1.5 × interquartile range (whiskers). (I) GSEA of autophagy of the HSF1-KD group 1, 2, and 12 h post-irradiation compared with the HSF1-OE group. (J) Correlation of HSF1 and ssGSEA scores of autophagy and cell cycle. (K) Differential expression of key proteins involved in these pathways 1, 2, and 12 h post-irradiation in the HSF1-OE/KD groups (left). Pearson’s correlation of these proteins with radiation dose (right). To elucidate the potential molecular mechanism by which HSF1 overexpression induced radioresistance, we performed proteomic comparative analysis between the HSF1-OE or HSF1-KD and control CRC cells in pre-irradiation and post-irradiation treatment with a series of irradiation doses and recovery times. Proteomic identifications showed no major differences in coverage between the HSF1-OE or HSF1-KD and control cells ([177]Figure S4D). We first compared proteome profiles of the HSF1-OE and control cells, and the results showed that 274 significantly up-regulated GPs in the HSF1-OE group were mainly enriched in VEGF, signaling by Rho GTPases, and cell cycle pathways (p < 0.05) ([178]Figure S4E; [179]Table S5). On the contrary, we found that these pro-proliferation-related signals were down-regulated in the HSF1-KD group ([180]Figure S4F). Further differential proteomic analysis between the HSF1-OE group/control group and the HSF1-KD group/control group validated that HSF1 overexpression promoted DNA replication and cell cycle, while HSF1 inhibition resulted in apoptosis and DNA damage response ([181]Figure 5D). The post-irradiation proteome data of CRC cells suggested that HSF1 overexpression rendered CRC cells prone to resistance to radiation (e.g., the upregulation of the cell cycle, DNA replication, and chromatin organization) (p < 0.05). In contrast, in the post-irradiation HSF1-KD CRC cells, we observed higher levels of autophagy- and lysosome-related signals, in addition to TNF alpha and apoptosis bioprocesses (p < 0.05). These results demonstrated that HSF1 overexpression increased resistance to radiation, and HSF1 inhibition could induce autophagy-related apoptosis ([182]Figure 5E; [183]Table S5). Next, we focused on the longitudinal trajectories of post-irradiation CRC cells and identified dynamic change of proteins during different recovery times in response to radiation. We compared the proteome data of post-irradiation HSF1-OE CRC cell within 1, 2, and 12 h recovery periods with pre-irradiation cells. We found that in the post-irradiation 1 h period, the main enriched bioprocesses were metabolism of lipids/nucleotides/RNA and vesicle-mediated transport, which provided energy and signaling molecules needed for proliferation.[184]^35 When cells proceeded into 2 h recovery after radiation, we observed the pathway enrichment of chromatin organization, suggesting that histone modification-related proteins were recruited to involve in chromatin organization and remodeling, which was reported to have an impact on the DNA-damage response and be modulated in response to DNA damage.[185]^36^,[186]^37^,[187]^38 In the post-irradiation 12 h period, the up-regulated proteins were enriched mainly in cell cycle regulation, DNA synthesis, and replication. These results revealed the dynamic signaling axis of metabolism of lipid-chromatin organization-cell cycle regulation occurred in HSF1-OE CRC cell, which was associated with radioresistance ([188]Figure 5F). Further ssGSEA revealed a gradual increase in cell cycle ssGSEA score (p < 1.0E-4, ANOVA) in the post-irradiation HSF1-OE group and a decrease in the post-irradiation HSF1-KD group. On the contrary, autophagy ssGSEA score showed a gradual decrease in the post-irradiation HSF1-OE group and an increase in the post-irradiation HSF1-KD group ([189]Figures 5G and 5H). Some studies showed that autophagy could enhance the anticancer effects of radiotherapy,[190]^39 which contributed to cell killing in response to radiation (defined as autophagic cell death).[191]^40 As shown in [192]Figure 5I, autophagy was definitely significantly enriched in the HSF1-KD group compared with the HSF1-OE group in the post-irradiation 1, 2, and 12 h periods. This finding validated HSF1 inhibition prone to resulting in autophagic cell death with radiation. In summary, HSF1 expression was negatively associated with autophagy ssGSEA score (Pearson’s r = −0.34, p = 3.1E-3) and positively associated with cell cycle ssGSEA score (Pearson’s r = 0.76, p = 6.9E-15) after radiation; a negative correlation was observed between cell cycle and autophagy ssGSEA score (Pearson’s r = −0.46, p = 4.7E-5) ([193]Figures 5J and [194]S4H). Specifically, we observed that autophagy markers (ATG7 and ATG9A) and apoptosis regulators (CASP6 and PARP2) were consistently down-regulated in the HSF1-OE group from 1 to 12 h post-irradiation and positively associated with irradiation dose, while lipid metabolism-related enzymes (AACS and DHCR7), chromatin organization-related enzymes (KDM2A and SETD2), and cell cycle proteins (CDCA5 and MCM4) were up-regulated in the HSF1-OE group in 1, 2, and 12 h post-irradiation, respectively ([195]Figure 5K). Overall, we propose that HSF1 overexpression drives a signaling axis of metabolism of lipid-chromatin organization-cell cycle regulation, thus triggering resistance to radiotherapy, while HSF1 inhibition is prone to produce autophagic cell death, thus promoting sensitivity to radiotherapy ([196]Figure 5L). Overexpression of hub protein HDAC6 regulated by NFYC synergistically induced cetuximab resistance With the increasing use of targeted therapies, more studies related to several trials have shown better clinical efficacy of cetuximab and bevacizumab compared with CT in KRAS wild-type and KRAS mutation metastatic CRC. However, the ORRs of cetuximab and bevacizumab were only about 50%. Therefore, it is imperative to explore the relative resistance/sensitivity mechanisms and identify the predictive markers for these targeted therapies. Here, we constructed a CTT subcohort (N = 31), in which patients receiving FOLFOX plus cetuximab therapy were most common (N = 17; N [S] = 9, N [NS] = 8). For patients treated with FOLFOX plus cetuximab therapy, we identified 64 proteins up-regulated in the S group and 47 proteins up-regulated in the NS group (p < 0.05, BH-adjusted p < 0.05) ([197]Figures 6A, [198]S5A, and S5B; [199]Table S6). Pathway enrichment revealed that the pathways PTEN-dependent cell-cycle arrest and apoptosis, oxidative phosphorylation, and neutrophil degranulation were significantly enriched in the S group, while cell-cycle-related pathways including the TNF alpha pathway, mRNA splicing, RNA transport, and DNA repair were significantly enriched in the NS group ([200]Figure S5C). Consistently, we found the reported EGRF-related signatures related to cetuximab sensitivity, such as PTEN[201]^41 and NF1,[202]^42 showed significant increases in the S group, while the reported signatures related to cetuximab resistance, mainly involving the downstream pathways of EGFR signaling, such as AKT2, MAPK13,[203]^43 MAP2K1,[204]^44 MTOR,[205]^45 and EPHA2,[206]^46 had relative increases in the NS group ([207]Figure S5D). Furthermore, we constructed protein-protein interaction (PPI) networks and identified hub proteins, of which BAG2 and HDAC6 were the core hub proteins dominant in the S and NS groups, respectively ([208]Figure 6B; [209]STAR Methods). Therefore, we proposed that BAG2 and HDAC6 were the important hub proteins associated with cetuximab sensitivity/resistance. Figure 6. [210]Figure 6 [211]Open in a new tab Differential expression of proteins and pathways in sensitive and non-sensitive groups of FOLFOX plus cetuximab subcohort (A) Differential proteins between the S and NS groups in the FOLFOX plus cetuximab subcohort. (B) PPI networks of hub proteins up-regulated in the S/NS groups. (C) Differential expression of HDAC6 in the S and NS groups of the subcohort in the Fudan cohort (S, N = 9; NS, N = 8) and [212]GSE59857 dataset (S, n = 10; NS, n = 9). (D) Diagram showing proteins with positive or negative correlation with HDAC6 and these proteins’ involved pathways. (E) GSEA enrichment of the proteins positively associated with HDAC6. (F) Differential expression of TFs in the S and NS groups. (G) Correlation of TFs with HDAC6. (H) Association of NFYC and HDAC6 on the basis of DepMap data. (I) Differential expression of HDAC6 and NFYC (S, N = 23; NS, N = 22) (left) and their correlation (right). (J) Differential expression of NFYC between the S and NS groups in the [213]GSE59857 dataset and the correlation of NFYC and HDAC6. (K) Qualification of HDAC6 stained by IHC in representative examples. The scale bar indicates 100 μm. Data are analyzed using two-sided Student’s t test and shown as mean ± SD (n = 3 independent experiments). (L) Schematic illustration indicating the potential regulation of HDAC6 in response to cetuximab treatment. To validate the association of HDAC6 with cetuximab resistance, we analyzed the transcriptional and pharmacological profiles of cell lines derived from CRC tissues from CRC patients treated with the same therapy regimen (FOLFOX plus cetuximab) from the [214]GSE59857 dataset ([215]STAR Methods).[216]^47 The results showed that the expression of HDAC6 was significantly increased in the NS group in the independent cohort (p = 0.017, Wilcoxon rank-sum test), consistent with the finding in our cohort ([217]Figure 6C). Then we performed correlation analysis and identified a group of proteins with positive correlations with HDAC6 (n = 4,817) (r > 0), which were involved mainly in proliferation-related pathways ([218]Figure 6D). GSEA also validated that the proteins positively associated with HDAC6 were significantly enriched in cell cycle (p = 0.033) ([219]Figure 6E). Matching with the Human Protein Atlas database, we found that 160 TFs were identified with positive correlations with HDAC6 ([220]Figure 6D). To explore the potential transcriptional regulation mechanism of cetuximab resistance related to HDAC6, we compared the proteome profiling of TFs. We concluded that 4 TFs were significantly differentially expressed between the S and NS groups, among which NFYC showed the highest correlation with HDAC6 ([221]Figures 6F and 6G). On the basis of CRC cell line data from DepMap, we observed consistent result (r = 0.51, p = 6.9E-6), which further validated the association between NFYC and HDAC6 ([222]Figure 6H). We also conducted PRM analysis in an independent cohort including 45 patients receiving FOLFOX plus cetuximab (N [S] = 23, N [NS] = 22). Consistently, NFYC and HDAC6 were also significantly elevated in the NS group (FC [NS/S] > 2; p < 0.05, Wilcoxon rank-sum test) and had a significantly positive association (r = 0.55, p = 1.2E-4) ([223]Figure 6I). In addition, we also validated that NFYC was up-regulated in the NS group in the [224]GSE59857 dataset ([225]STAR Methods), and correlation analysis demonstrated a significant positive association between NFYC and HDAC6 (r = 0.48, p = 0.04) ([226]Figure 6J). Furthermore, we performed analysis by the deposited matrix of binding domain for NFYC in Ensembl on the basis of the JASPAR open-access database ([227]https://jaspar.genereg.net). JASPAR analysis revealed that HDAC6 carried 2 putative NFYC-binding sites (GGCCAATAATA and TTCCACTCAGA) along its DNA transcriptional regulatory region, with homology higher than 80% ([228]Table S6), which provided the potential indirect transcriptional inhibition target for cetuximab targeted therapy. More directly, IHC showed significant differences in HDAC6 expression between the S and NS groups in patients treated with cetuximab (p < 0.0001, Student’s t test) ([229]Figure 6K), demonstrating the direct inhibition target for cetuximab-resistant patients. Taken together, HDAC6 regulated by the TF NFYC could synergistically resulted in cetuximab resistance by inducing tumor proliferation ([230]Figure 6L). Therefore, we propose a potential strategy whereby cetuximab plus HDAC6 inhibitors (such as tubastatin A, ACY-241, and nexturastat), which could improve the response of cetuximab therapy in the clinic. Construction and validation of predictive models for CRC CT, chemoradiation, and targeted therapies Having proposed the potential sensitive or resistant mechanism for the personalized therapy guidance, we next set out to determine whether DEPs could distinguish S patients from NS patients in response to FOLFOX therapy, FOLFOX plus radiotherapy, and FOLFOX plus cetuximab ([231]Figure 7A). On the basis of these DEPs, we identified a subset of signatures that accurately discriminates S and NS patients in response to these therapies by stepwise logistic regression. To train and subsequently test the model, samples were partitioned on the basis of response group, and 80% of samples were used as a training set, with the remaining 20% representing the independent testing set. On the basis of these S/NS signatures (FOLFOX, n = 5 [LPCAT3, UNC13B, VWA1, NDUFB3, and SIRT5]; FOLFOX plus radiotherapy, n = 5 [MOXD1, PFN2, DERL2, SHMT2, and MTHFD2]; FOLFOX plus cetuximab, n = 4 [CES2, PRKG1, NUP85, and HDAC6]), we applied 10-fold cross-validation to the training set to evaluate the performance ([232]STAR Methods). Further accuracy metrics demonstrated that the three predictive models had high accuracy of 0.97, 0.88, and 0.85 in the FOLFOX therapy, FOLFOX plus radiotherapy, and cetuximab therapy subcohorts, respectively ([233]Figure 7B). When applied to the independent test set samples, the three predictive models separately achieved high accuracy of 0.88, 1, and 0.75 in the three subcohorts ([234]Figure 7B). For further testing the predictive power of three models for therapy response, we then applied the predictive models for FOLFOX and FOLFOX plus radiotherapy in the CT and CRT subcohorts, which achieved high accuracy of at least 0.82 ([235]Figure S5E). Furthermore, the predictive models for FOLFOX and FOLFOX plus cetuximab therapies were validated in other public datasets, achieving high accuracy of 0.76, 0.93, and 1 in the [236]GSE19860 dataset, the [237]GSE5851 dataset, and Project ID [238]OEP000729 ([239]STAR Methods), respectively ([240]Figure S5E). In addition, to verify the stability of the predictive performance, we also applied 3-fold cross-validation to the training sets of three models for the FOLFOX, FOLFOX plus radiotherapy, and FOLFOX plus cetuximab therapy subcohorts and then compared the predictive effect of 10-fold cross-validation and 3-fold cross-validation ([241]STAR Methods). As a result, the average AUC values of predictive models with 10-fold cross-validation and 3-fold cross-validation on training sets were 0.90 and 0.81 (FOLFOX), 0.92 and 0.82 (FOLFOX plus radiotherapy), and 0.81 and 0.77 (FOLFOX plus cetuximab therapy), which showed a comparable predictive effect between 3-fold and 10-fold cross-validation among the predictive models for FOLFOX, FOLFOX plus radiotherapy, and FOLFOX plus cetuximab ([242]Figure S5F). Overall, the predictive models in the discovery cohort exhibited robustness, accuracy, and stability in the external validation cohort. Figure 7. [243]Figure 7 [244]Open in a new tab Construction and validation of predictive classifiers (A) Diagram describing construction and validation of predictive modelling for CRC therapy. (B) Heatmap showing differential expression of predictive model signatures. Classification error matrix of 80% training set and 20% testing set in the three subcohorts. (C) Heatmaps, classification error matrix, and ROC curves showing the differential expression of signatures quantified using the PRM approach in three independent subcohorts (two-sided Wilcoxon rank-sum test; ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001). To evaluate the accuracy of the predictive signatures for the therapeutic response, we designed a PRM strategy to quantify these S/NS signatures in the independent cohort composed of 130 CRC patients receiving FOLFOX (N = 45; N [S] = 18, N [NS] = 27), FOLFOX plus radiotherapy (N = 40; N [S] = 21, N [NS] = 19), or FOLFOX plus cetuximab (N = 45; N [S] = 23, N [NS] = 22) therapies ([245]Table S7). On the basis of the PRM quantification of these S/NS signatures ([246]STAR Methods), we performed comparative analysis of S/NS signatures between the S and NS groups in three independent subcohorts, respectively. As a result, we observed the significantly differential expression of these signature proteins between the S and NS groups in three independent subcohorts (FC > 2; p < 0.05, Wilcoxon rank-sum test) ([247]Figure S5G; [248]Table S7). The heatmaps showed a clear separation between the S and NS groups in the three independent subcohorts ([249]Figure 7C). In addition, we further validated the predictive models by accuracy metrics and receiver operating characteristic (ROC) analysis on the basis of the PRM quantification data, achieving high accuracy of 0.93, 0.93, and 0.96 in the FOLFOX, FOLFOX plus radiotherapy, and FOLFOX plus cetuximab subcohorts, respectively. Moreover, the ROC curves show high sensitivity and specificity of prediction, with AUC values of 0.99, 0.94, and 0.99 in these subcohorts, respectively ([250]Figure 7C). Collectively, the predictive power of these signature proteins was further validated in the independent validation cohort by PRM assays. In addition, these predictive models exhibited good performance in CT, CRT, and CTT, demonstrating robustness, accuracy, and universality. Discussion According to the estimated number of total deaths, the mortality rate of CRC was ranked fourth in China.[251]^48 CT/CRT has been advocated as the standard treatment for CRC.[252]^49 With the development of targeted drug agents, such as VEGF inhibitors and EGFR inhibitors, combination with chemotherapies could improve overall survival and have been approved for clinical use. There are no indicators to predict these therapies’ effectiveness, and the development of drug resistance remains unresolved. Therefore, it is imperative to make a preclinical diagnosis of tumor response to these first-line therapies. Herein, we present a large-scale clinical proteomic landscape, aiming to identify reliable predictive markers of the CRC patient response to the diverse therapies. In this study, proteomic subtypes were identified with distinct molecular features and therapy response. Pathway enrichment analysis suggested immune activation or RAS-MEK-ERK inhibition could be applied in the combination of CT and targeted therapy, thus improving therapy sensitivity and reducing therapy resistance. The key findings derived from proteomic subtypes deepened our understanding of the landscape, including different therapy regimens, providing the opportunity to link our subtypes with other molecular subtypes identified in other public independent cohorts and important clues for us to validate in the further analysis. Further comparative analysis and exploration between S patients and NS patients provided extensive resistance/sensitivity mechanisms related to different therapies. To further explore the biological alterations associated with therapy response in the post-therapy period, we collected 30 post-therapy samples matched with their paired pre-therapy samples from patients receiving FOLFOX (N [S] = 5, N [NS] = 5), FOLFOX plus radiotherapy (N [S] = 5, N [NS] = 5), or FOLFOX plus cetuximab therapy (N [S] = 5, N [NS] = 5). GSEA results revealed the dominant pathways correlated with therapy response. As a result, we found that the S group of patients after receiving FOLFOX therapy was enriched by IL-2 signaling pathway, MHC pathway, interferon gamma signaling, and autophagy, while the NS group was dominant with proliferation-related pathways and glycerophospholipid/pyrimidine metabolism, which showed similar enrichment of the S and NS groups in the post-therapy and pre-therapy samples originated from CRC patients receiving FOLFOX therapy ([253]Figure S6A; [254]Table S8). In addition, we also found the similar enrichment associated with therapy response in the pre-therapy and post-therapy samples from patients receiving FOLFOX plus radiotherapy or cetuximab therapy ([255]Figures S6B and S6C; [256]Table S8). Moreover, we also observed that the key proteins more highly expressed in the pre-therapy NS group compared with the S group exhibited an enhanced high-expression tendency in the post-therapy NS group. In contrast, the key proteins more highly expressed in the pre-therapy S group compared with the NS group exhibited a weakened high-expression tendency in the post-therapy S group. Taken together, these results suggested an enhancement of the negative association of molecular features with therapy response but an attenuation of the positive association in the post-therapy samples compared with the pre-therapy samples ([257]Figure S6D; [258]Table S8). In this study, we integrated pre-therapy samples with response evaluation to explore drug sensitivity and resistance, which could well reveal the original biological heterogeneity in the tumor themselves of patients without being affected by the secondary effect of treatment. More important, we developed predictive models with high accuracy for each therapy in the discovery cohort and validated the prediction performance in the independent validation cohort through PRM quantification. In this study, we should address that these tissue samples originated from therapy-naive patients, which indicated that the original heterogeneity in the tumor themselves of different patients; such biomarkers could be used for prediction analysis to determine which treatment is effective for the patients in the future and to stratify patients into specific subcategories of therapy response, thus avoiding the use of inappropriate treatment and missing the best time for treatment, ultimately improving the clinical effectiveness and the quality of clinical strategy design. Collectively, this study provides potential molecular features associated with therapeutic response, implicates the underlying regulatory mechanisms, and identifies biomarkers for therapy response, ultimately serving to benefit CRC patients in clinical practice. Finally, we integrated key findings to present a comprehensive view ([259]Figure S6E), revealing response heterogeneity for different therapies and improving the therapeutic effectiveness and development of personalized therapy for CRC patients. The construction and validation of proteomic-based predictive classifiers for these therapies of CRC contributed to personalized therapies to a certain extent, thereby moving us toward the era of proteomic-driven precision medicine. More important, we propose that inhibition of HSF1 can facilitate more effective clinical decision making to determine the relative benefit of radiotherapy for patients. In summary, our study has potentially important clinical implications in CRC, and we look forward to further developing potential therapeutic combinations for these therapies. Limitations of the study First, as a single-center study, this work likely does not represent all the heterogeneous mechanisms underlying clinical resistance to first-line therapies. Second, the biomarkers identified in this study are based on proteome data, and the application of these biomarkers to clinical decision making is still to be validated. Third, this study focuses on the initial first-line therapy, and longitudinal follow-up or prospective studies are also needed. Fourth, clinical assays to measure these markers reliably in a routine method are yet to be developed. STAR★Methods Key resources table REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies __________________________________________________________________ Rabbit Polyclonal anti-HSF1 antibody Signalway Antibody Cat# 41037-1; RRID: [260]AB_3073724 Rabbit polyclonal antibody against HDAC6 Signalway Antibody Cat# 40999-1; RRID: [261]AB_3073725 anti-γ-H2AX antibody ABclonal Cat# AP0687; RRID: [262]AB_2863808 FITC-labeled Goat Anti-Rabbit IgG Beyotime Cat# A0562; RRID: [263]AB_2923335 Anti-β-Actin antibody abmart Cat# P30002M; RRID: [264]AB_2936505 __________________________________________________________________ Biological samples __________________________________________________________________ Archival formalin-fixed paraffin-embedded (FFPE) colorectal cancer samples Zhongshan Hospital, Fudan University This paper __________________________________________________________________ Chemicals, peptides, and recombinant proteins __________________________________________________________________ HPLC-grade water J.T. Baker Cat# 4218-03 Urea Sigma Cat# U5378 Sodium chloride Sigma Cat# S5886 Tris Amresco Cat# 0497 Phenylmethylsulfonyl fluoride Amresco Cat# M145 Dithiothreitol Sigma Cat# 43815 Formic acid Sigma Cat# 0507 Trypsin Promega Cat# V528A Ammonium bicarbonate J.T. Baker Cat# 3003-01 Acetonitrile J.T. Baker Cat# 9829-03 Hexadimethrine Bromide Sigma Cat# H9268 NONIDET P-40 SUBSTITUTE Amresco Cat# M158 Colchicine Sangon Biotech Cat# A600322-0100 Bradford Dye Reagent TaKaRa Cat# T9310A-1 SuperSignal West Pico PLUS Chemiluminescent Substrate Thermo Scientific Cat# 34580 Lipofectamine 2000 Invitrogen Cat# 11668-019 Colchicine Sangon Biotech Cat# A600322-0100 DAPI Sigma Cat# D8417 Q Exactive HF-X Hybrid Quadrupole-Orbitrap Mass Spectrometer ThermoFisher Cat# 0726042 EASY-nLC 1200 system ThermoFisher Cat# LC140 Amersham Imager 600 GE Healthcare Bio-Sciences AB Cat# 66720240 iMark Microplate Reader BioRad Cat# 15001 Vacuum centrifuge Eppendorf Cat# 07-748-15 __________________________________________________________________ Critical commercial assays __________________________________________________________________ Easypure Plasmid Miniprep Kit TransGen Biotech Cat# EM101-02 RevertAid First Strand cDNA Synthesis Kit Thermo Scientific Cat# K1622 Cell Counting Kit-8 Beyotime Cat# C0046 __________________________________________________________________ Deposited data __________________________________________________________________ Proteomic data of the colorectal cancer cohort This paper iProx: [265]IPX0002629000 ProteomeXchange: [266]PXD046566 The TCGA cohort Zhang et al.[267]^11 [268]https://cptac-data-portal.georgetown.edu Transcriptomic expression data of the rectal cancer patients Millino et al.[269]^32 GEO accession - [270]GSE68204 Global gene expression profiling of human colorectal cancer cell lines Medico et al.[271]^47 GEO accession - [272]GSE59857 Transcriptional data from CRC patients receiving cetuximab therapy Khambata-Ford et al.[273]^50 GEO accession - [274]GSE5851 Proteomic data from CRC patients receiving cetuximab therapy Li et al.[275]^13 [276]OEP000729 Gene expression profiles from CRC patients receiving FOLFOX therapy Gene Expression Omnibus GEO accession - [277]GSE19860 __________________________________________________________________ Experimental models: Cell lines __________________________________________________________________ Human: SW480 cell line ATCC RRID: CVCL_0546 Human: SW620 cell line ATCC RRID: CVCL_0547 Human: RKO cell line ATCC RRID: CVCL_0504 Human: HCT116 cell line ATCC RRID: CVCL_0291 __________________________________________________________________ Recombinant DNA __________________________________________________________________ plvx-PURO vector Takara Bio Cat# 632164 pMD2.G plasmid Addgene Cat# 12259 psPAX2 plasmid Addgene Cat# 12260 pLKO.1-puro vector Addgene Cat# 8453 __________________________________________________________________ Software and algorithms __________________________________________________________________ R version v 3.5.1 The R Project for Statistical Computing [278]https://www.r-project.org/ RStudio (v1.2.5019) R Studio Team, 2015 [279]https://rstudio.com/products/rstudio/download/ ConsensusClusterPlus (version 3.8) Monti; Wilkerson[280]^51^,[281]^52 [282]https://bioconductor.org/packages/release/bioc/html/ConsensusClust erPlus.html Survminer R package (version 0.2.4) NA [283]https://cran.r-project.org/web/packages/survminer/index.html Complex Heatmap R package (version 2.0.0) Gu[284]^53 [285]http://www.bioconductor.org/packages/release/bioc/html/ComplexHeat map.html pROC R package (version 1.16.2) NA [286]https://cran.r-project.org/web/packages/pROC/index.html Caret R package (version 6.0–86) NA [287]https://cran.r-project.org/web/packages/caret/index.html ImageJ software (version 1.52a) National Institutes of Health [288]https://imagej.nih.gov/ij/ GraphPad Prism (version 8.0.1) GraphPad [289]https://www.graphpad.com/scientific-software/prism/ Firmiana Feng[290]^54 [291]https://phenomics.fudan.edu.cn/firmiana/gardener/ Mascot (version 2.3) Matrix Science [292]http://www.matrixscience.com/ NCBI human Refseq protein database NCBI [293]https://www.ncbi.nlm.nih.gov/refseq/ RapidMiner 9.6.0 RapidMiner, Inc [294]https://rapidminer.com/ [295]Open in a new tab Resource availability Lead contact Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Chen Ding (chend@fudan.edu.cn). Materials availability This study did not generate new unique reagents. Data and code availability * • The raw mass spectrometry proteomics data generated in this study have been deposited in the ProteomeXchange Consortium (dataset identifier: [296]PXD046566) via the iProX partner repository ([297]https://www.iprox.cn/)[298]^55 under project ID [299]IPX0002629000. The public dataset related to CRC patients and CRC cell lines were downloaded from Gene Expression Omnibus (GEO, [300]https://www.ncbi.nlm.nih.gov/geo/), including the transcriptomic expression data from rectal cancer patients in the [301]GSE68204 dataset [[302]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE68204], the transcriptional and pharmacological profiles of stable cell lines derived from colorectal cancer and related tissues in the [303]GSE59857 dataset [[304]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE59857], the gene expression profiles from CRC patients receiving FOLFOX therapy in the [305]GSE19860 dataset [[306]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE19860], the transcriptional data from CRC patients receiving cetuximab therapy in the [307]GSE5851 dataset [[308]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE5851], and the proteomic data from CRC patients receiving cetuximab therapy in the Project ID [309]OEP000729 [[310]https://www.biosino.org/node/project/detail/OEP000729]. Source data are provided with this paper. * • This paper does not report original code. * • Any additional information required to reanalyze the data reported in this work paper is available from the lead contact upon request. Experimental model and study participant details Clinical sample acquisition This study was approved by the Research Ethics Committee of Zhongshan Hospital (B2019-200R). Written informed consent was received from all patients included in this study. Archival formalin-fixed, paraffin-embedded (FFPE) tissues from patients with colorectal cancer (CRC), from 2009 to 2021, were reviewed in the Department of Pathology, Zhongshan Hospital, Fudan University (Shanghai, R. P. China). The inclusion criteria for clinical samples in this study were as follows: (i) the tissue source was from untreated patients undergoing biopsy surgery for a diagnosis of CRC; (ii) areas containing 80% or more tumor were examined independently by three expert gastrointestinal pathologists (D.J., C.X. and Y.H); (iii) no history of severe heart or liver disease, psychiatric disorders, hemorrhage, or perforation of the digestive tract; (iv) including CRC patients receiving the first-line therapy (including chemotherapy, or chemoradiotherapy, or targeted therapy), excluding the CRC patients treated with the mixed treatments (such as second-line or even third-line therapy); (v) the complete therapy response evaluation information followed by the first therapy. According the strict inclusion criteria for CRC patient in this study, we finally collected 254 biopsy tumor FFPE samples derived from therapy-naïve CRC patients, including discovery cohort (N = 124) and validation cohort (N = 130), receiving either chemotherapy (CT, mainly FOLFOX therapy, composed of oxaliplatin, fluorouracil, and leucovorin), chemoradiotherapy (CRT, mainly FOLFOX plus radiotherapy), or chemotherapy combined with targeted therapy (CTT, mainly FOLFOX plus cetuximab). In the discovery cohort, the 124 patients were receiving chemotherapy (N = 48), chemoradiotherapy (N = 45), or chemotherapy combined with targeted therapy (N = 31). The independent validation cohort for PRM verification included 130 patients with CRC, receiving either FOLFOX therapy (N = 45), FOLFOX plus radiotherapy (N = 40), or FOLFOX plus cetuximab therapy (N = 45). In this study, all the preoperative therapies were the first-line therapies, and all the therapy regimens were given at standard dosing as described in previous studies.[311]^56^,[312]^57 The archival formalin-fixed paraffin-embedded (FFPE) tissues were collected from the therapy-naïve patients for proteome measurement (as known as pre-therapy samples), and histologically scored by three expert gastrointestinal pathologists (D.J., C.X. and Y.H). Medical records of each patient were reviewed to obtain the follow-up data, especially the therapy response. In this study, the treatment response was evaluated by CT/MRI scanning following the Response Evaluation Criteria in Solid Tumors (RECIST) (version1.1) guideline.[313]^58 Tumor response was assessed and categorized as a complete response (CR), partial response (PR), stable disease (SD), or progressive disease (PD). All the process was followed by this guideline, and we adopted the consistent imaging approach within each patient (before and after the treatment) and consistent evaluation standard across all patients during the therapy period. In clinic, the ORR, generally defined by the Food and Drug Administration as the sum of PR plus CR, is a direct measurement of drug antitumor activity. Here, the ORR, defined as PR plus CR, was selected for the efficacy evaluation; patients with CR and PR were defined as sensitive (S) and those with SD and PD were defined as non-sensitive (NS). Consistently, the distinction and definition has been adopted in the similar published researches.[314]^59^,[315]^60 The chemotherapy (CT) subcohort was grouped into CT-sensitive (S) group (CT (S), N = 19) and CT-non-sensitive (NS) group (CT (NS), N = 29), of which, FOLFOX subcohort was composed of 16 S and 22 NS; the chemoradiotherapy (CRT) subcohort was grouped into CRT-sensitive group (CRT (S), N = 27) and CRT-non-sensitive group (CRT (NS), N = 18), of which, FOLFOX + radiotherapy subcohort was composed of 19 S and 11 NS; the chemotherapy combined with targeted therapy (CTT) subcohort was grouped into CTT-sensitive group (CTT (S), N = 18) and CTT-non-sensitive group (CTT (NS), N = 13), of which, FOLFOX + cetuximab subcohort was composed of 9 S and 8 NS. The independent validation cohort for PRM verification composed with patients receiving FOLFOX therapy (N = 45, N (S) = 18; N (NS) = 27), patients receiving FOLFOX + radiotherapy (N = 40, N (S) = 21; N (NS) = 19), patients receiving FOLFOX + cetuximab therapy (N = 45, N (S) = 23; N (NS) = 22). All the FFPE tissues with at least 80% tumor purity were collected from a total of 254 therapy-naïve patients with CRC. The detailed clinical characteristics and definitions were shown in [316]Table S1; [317]Figures S1A and S1B. Cell line Human HEK293T (Cat# CRL-3216 from ATCC; RRID: CVCL_0063), RKO (Cat# CRL-2577 from ATCC; RRID: CVCL_0504), HCT116 (Cat# CRL-247 from ATCC; RRID: CVCL_0291), SW480 (Cat# CCL-228 from ATCC; RRID: CVCL_0546), and SW620 (Cat# CCL-227 from ATCC; RRID: CVCL_0547), were obtained and cultured in DMEM (GIBCO) with 10% FBS (GIBCO) in 5% CO[2] at 37°C. Cell lines were examined by short tandem repeat markers (STR) profiling and then analyzed using the DSMZ (German Collection of Microorganisms and Cell Cultures) online STR database ([318]http://www.dsmz.de/fp/cgi-bin/str.html). Cell lines were tested negative for mycoplasma contamination and grown according to the instruction. Method details Protein extraction and trypsin digestion of CRC FFPE samples The biopsy tumor FFPE samples derived from therapy-naïve CRC patients were collected, and the tumor regions were determined by pathological examination. For clinical sample preparation, sections with 10 μm thick and 10 mm^2 area from FFPE blocks were macro-dissected, deparaffinized with xylene, and washed with ethanol. The ethanol was removed completely and the sections were left to air-dry. For this purpose, a hematoxylin-stained section of the same tumor was used as reference. Areas containing 80% or more tumor were examined independently by three expert gastrointestinal pathologists (D.J., C.X. and Y.H). Lysis buffer [0.1 M Tris-HCl (pH 8.0), 0.1 M DTT (Sigma, 43815), 1 mM PMSF (Amresco, M145)] was added to the same weight of tissues (around 3,470 μg), and subsequently sonicated for 1 min (3 s on and 3 s off, amplitude 25%) on ice. The supernatants were collected, and the protein concentration was determined using the Bradford assay. The extracted tissues were then lysed with 4% sodium dodecyl sulfate (SDS) and kept for 2–2.5 h at 99°C with shaking at 1,800 rpm. The solution was collected by centrifugation at 12,000 ×g for 5 min. A 4-fold volume of acetone was added to the supernatant and kept in −20°C for a minimum of 4 h. Subsequently, the acetone-precipitated proteins were washed three times with cooled acetone. Filter-aided sample preparation (FASP) procedure was used for protein digestion.[319]^61 The proteins were resuspended in 200 μL 8 M urea (pH 8.0) and loaded in 30 kD Microcon filter tubes (Sartorius) and centrifuged at 12,800 g for 20 min. The precipitate in the filter was washed three times by adding 200 μL 50 mM NH[4]HCO[3]. The precipitate was resuspended in 50 μL 50 mM NH[4]HCO[3]. Protein samples underwent trypsin digestion (enzyme-to-substrate ratio of 1:50 at 37°C for 18–20 h) in the filter, and then were collected by centrifugation at 12,800 g for 15 min. Additional washing, twice with 200 μL of water, was essential to obtain greater yields. In the sample preparation of this study, around 320 μg proteins (9.25% yield) were extracted for around 3,470 μg tissues of each sample with a linear relationship ([320]Figures S1D and S1E), and ultimately a standard loading amount of 500 ng peptide for the proteomic measurement. Finally, the centrifugate was pumped out using the AQ model Vacuum concentrator (Eppendorf, Germany). LC-MS/MS Peptide samples were analyzed on a Q Exactive HF-X Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Fisher Scientific, Rockford, IL, USA) coupled with a high-performance liquid chromatography system (EASY-nLC 1200, Thermo Fisher Scientific). Peptides, re-dissolved in Solvent A (0.1% formic acid in water), were loaded onto a 2-cm self-packed trap column (100-μm inner diameter, 3-μm ReproSil-Pur C18-AQ beads, Dr. Maisch GmbH) using Solvent A, and separated on a 150-μm-inner-diameter column with a length of 15 cm (1.9-μm ReproSil-Pur C18-AQ beads, Dr. Maisch GmbH) over a 75 min gradient (Solvent A: 0.1% formic acid in water; Solvent B: 0.1% formic acid in 80% ACN) at a constant flow rate of 600 nL/min (0–75 min, 0 min, 4% B; 0–10 min, 4–15% B; 10–60 min, 15–30% B; 60–69 min, 30–50% B; 69–70 min, 50–100% B; 70–75 min, 100% B). The eluted peptides were ionized under 2 kV and introduced into the mass spectrometer. MS was operated under a data-dependent acquisition (DDA) mode. For the MS1 Spectra full scan, ions with m/z ranging from 300 to 1,400 were acquired by Orbitrap mass analyzer at a high resolution of 120,000. The automatic gain control (AGC) target value was set as 3E+06. The maximal ion injection time was 80 ms. MS2 Spectra acquisition was performed in top-speed mode. Precursor ions were selected and fragmented with higher energy collision dissociation with a normalized collision energy of 27%. Fragment ions were analyzed using an ion trap mass analyzer with an AGC target value of 5E+04, with a maximal ion injection time of 20 ms. Peptides that triggered MS/MS scans were dynamically excluded from further MS/MS scans for 12 s. A single-run measurement was kept for 75 min. All data were acquired using Xcalibur software v2.2 (Thermo Fisher Scientific). Peptide and protein identification MS raw files were processed using the Firmiana proteomics workstation.[321]^54 Briefly, raw files were searched against the NCBI human Refseq protein database (released on 04-07-2013; 32,015 entries) using the Mascot search engine (version 2.4, Matrix Science Inc). The mass tolerances were: 20 ppm for precursor and 50 mmu for product ions collected by Q Exactive HF-X. Up to two missed cleavages were allowed. The database searching considered cysteine carbamidomethylation as a fixed modification, and N-acetylation, and oxidation of methionine as variable modifications. Precursor ion score charges were limited to +2, +3, and +4. For the quality control of protein identification, the target-decoy-based strategy was applied to confirm the FDR of both peptide and protein, which was lower than 1%. Percolator was used to obtain the quality value (q-value), validating the FDR (measured by the decoy hits) of every peptide-spectrum match (PSM), which was lower than 1%. Subsequently, all the peptides shorter than seven amino acids were removed. The cutoff ion score for peptide identification was 20. All the PSMs in all fractions were combined to comply with a stringent protein quality control strategy. We employed the parsimony principle and dynamically increased the q-values of both target and decoy peptide sequences until the corresponding protein FDR was less than 1%. Finally, to reduce the false positive rate, the proteins with at least one unique peptide were selected for further investigation. Label-free-based MS quantification of proteins The one-stop proteomic cloud platform “Firmiana” was further employed for protein quantification. Identification results and the raw data from the mzXML file were loaded. Then for each identified peptide, the extracted-ion chromatogram (XIC) was extracted by searching against the MS1 based on its identification information, and the abundance was estimated by calculating the area under the extracted XIC curve. For protein abundance calculation, the nonredundant peptide list was used to assemble proteins following the parsimony principle. The protein abundance was estimated using a traditional label-free, intensity-based absolute quantification (iBAQ) algorithm,[322]^62 which divided the protein abundance (derived from identified peptides’ intensities) by the number of theoretically observable peptides. A match between runs was enabled to transfer the identification between separate LC-MS/MS runs based on their accurate mass and retention time after retention time alignment.[323]^63 We built a dynamic regression function based on the commonly identified peptides in tumor samples. According to correlation value R^2, Firmiana chose linear or quadratic functions for regression to calculate the retention time (RT) of corresponding hidden peptides, and to check the existence of the XIC based on the m/z and calculated RT. Subsequently, the fraction of total (FOT), a relative quantification value was defined as a protein’s iBAQ divided by the total iBAQ of all identified proteins in one experiment, and was calculated as the normalized abundance of a particular protein among experiments. Finally, the FOT values were further multiplied by 10^5 for ease of presentation, and missing values were assigned 10^−5.[324]^15 Quality control of the mass spectrometry data To quality control the MS performance, the HEK293T cell lysate was measured every 20 samples as the quality control standard, which was adopted in proteomic studies.[325]^15^,[326]^16^,[327]^64 The quality control standard was digested and analyzed using the same method and conditions as the CRC samples. Pearson’s correlation coefficient was calculated for all quality control runs using the R statistical analysis software v.3.5.1 ([328]Figure S1C). The average correlation coefficient among the standards was 0.971, and the maximum and minimum values were 0.98 and 0.95, respectively. These results demonstrated the consistent stability of the MS platform. The log[10] transformed FOTs for each CRC sample ([329]Figure S1F) were plotted to show consistency of data quality. Targeted PRM analysis Using the library search results, a set of target peptides that unique to proteins was selected and parallel reaction monitoring (PRM) method was designed, including immune-regulated proteins (including RNASE2, PTPRJ, CEACAM8, SERPINA5, C7, and F10), and S/NS signatures (LPCAT3, UNC13B, VWA1, NDUFB3 and SIRT5 included in the predictive model for FOLFOX, MOXD1, PFN2, DERL2, SHMT2 and MTHFD2 included in the predictive model for FOLFOX + radiotherapy, and CES2, PRKG1, NUP85 and HDAC6 included in the predictive model for FOLFOX + cetuximab therapy). Besides, house-keeping proteins, such as VCP, RPLP0, PSMB4, were also included for the reference. Equal amount of tumor tissue from each sample (an independent cohort composed of 130 CRC patients, including 62 sensitive patients and 68 non-sensitive patients) was digested as described in the part of profiling preparation. Peptide samples were injected into the Q Exactive HF-X Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Scientific) operating in PRM mode with quadrupole isolation and HCD fragmentation. The full MS mode was measured at resolution 60,000 with AGC target value of 3E6 and maximum IT of 20 ms, with scanning range of 300–1400 m/z. Target ions were submitted to MS/MS in the HCD cell (1.6m/z isolation width, 27% normalized collision energy). 130 PRM events were performed after MS1 scanning, at resolution 15,000 with AGC target value of 1E6 and maximum IT of 25 ms. Separation was achieved on a 150-μm-inner-diameter column with a length of 15 cm (1.9-μm ReproSil-Pur C18-AQ beads, Dr. Maisch GmbH) in an Easy 1200 nLC HPLC system (Thermo Scientific). Solvent A was 0.1 formic acid in water and solvent B was 0.1% formic acid, 80% ACN in water. Peptides were separated at 600 nL/min across a gradient ranging from 4% to 100% B over 75 min (0–75 min, 0 min, 4% B; 0–10 min, 4–15% B; 10–60 min, 15–30% B; 60–69 min, 30–50% B; 69–70 min, 50–100% B; 70–75 min, 100% B). Raw data was searched by Skyline-daily (4.2.1.19004, University of Washington, USA). The proteins were quantified with the fragment total area reported by Skyline-daily. We selected peptides and tested their stability of signal and shape of peaks in the pool sample for final quantification, and referred to the ranking offered by skyline. Construction of HSF1 overexpression, lentivirus production and cell transduction For the analysis of HSF1 in tumor proliferation and tumor response to radiation, we constructed stable cell lines overexpressing HSF1-FLAG in two CRC cells with low HSF1 expression (SW480 and SW620). The cDNA of HSF1 was cloned into the plvx-PURO vector via the unique XhoI site and the neighboring EcoRI site. For convenient detection, an FLAG tag encoding sequence (GATTACAAGGATGACGACGATAAG) was inserted before the start codon (ATG) to express the HSF1-FLAG fusion protein. The primers used for plasmid construction as following: Forward primer (5′–3′): ctcgagatggattacaaggatgacgacgataagATGGATCTGCCCGTGGGCCCC, Reverse primer (5′–3′): GAATTCCTAGGAGACAGTGGGGTCCTT. The double-stranded HSF1 DNA was cloned into the plvx-PURO vector; 4 μg of packaging plasmids psPAX2: pMD2.G (3:1) and 4 μg of Lenti-vector containing the target gene were co-transfected into 2.5×10^6 HEK-293T cells using Polyethylenimine (PEI). In this transfection method, 500 μL of DMEM (serum-free medium) and the plasmid were placed in an empty EP tube and PEI (three times the concentration of plasmid) was added into the medium with vigorous shaking. The mixture was incubated for 15 min. Meanwhile, the cell culture medium was replaced with 2 mL of fresh 10% FBS medium. After 15 min, the mixture was added to the cells, and the fresh medium was replaced after 8 h. After 48–72 h, the transfection was completed and the cells were treated or harvested. The media containing the lentivirus particles were collected after 48–72 h, and centrifuged at 5,000 xg for 5 min. The supernatants approximately 1.5 mL mixed with 0.5 mL fresh medium were used independently to infect SW480 and SW620 cells; 24 h later, removed the media containing the lentivirus particles and added the fresh medium in the presence of hexadimethrine bromide (Sigma, H9268) at final concentration of 1 μg/mL for 12 h. After infection, cells were cultured and selected with puromycin for the generation of stably overexpressed cells. Empty plvx-PURO vector was used as a negative control. Construction of HSF1 shRNA, lentivirus production and cell transduction We also constructed HSF1-knockdown cell line (HSF1-KD) in two CRC cells with high HSF1 expression (HCT116 and RKO). The pLKO.1-puro vectors carrying one of two different HSF1 shRNA constructs were custom-designed by using online software available from the Sigma-Aldrich. HSF1 shRNA sequences were 5′-CCAGCAACAGAAAGTCGTCAA-3′. A nonsense sequence shRNA (5′-CAACAAGATGAAGAGCACCAA-3′) was used as a negative control. All these shRNA sequences were prepared and inserted into the pLKO.1-puro vector digested by EcoRI and AgeI. All plasmids were synthesized and amplified and confirmed by sequencing prior to use. Correct insertions of shRNA cassettes were confirmed by restriction mapping and direct DNA sequencing. The shRNA-expressing lentiviral were transfected into 293T cells using Polyethylenimine (PEI). Infectious lentiviruses were harvested 48–72 h post-transfection, centrifuged to remove cell debris. The supernatants approximately 1.5 mL mixed with 0.5 mL fresh medium were used independently to infect RKO and HCT116 cells; 24 h later, removed the media containing the lentivirus particles and added the fresh medium in the presence of hexadimethrine bromide (Sigma, H9268) at final concentration of 1 μg/mL for 12 h. After infection, cells were cultured and selected with puromycin for the generation of stably overexpressed cells. Empty pLKO.1-puro vector was used as a negative control. The result related to HSF1 shRNA2 was shown in the manuscript. Proteome profiling and differential analysis of CRC cell lines The protein concentration was determined using the Bradford assay. Cells were boiled in a 99°C metal bath for 30 min with 100 μL of 50 mM ABC buffer (ammonium bicarbonate) containing SDS at a final concentration of 4%. Protein samples underwent trypsin digestion (enzyme-to-substrate ratio of 1:50 at 37°C for 18–20 h). All peptide samples were desalinated using a C18 column (50% acetonitrile and 0.1% formic acid) and then analyzed using a Q Exactive HF-X Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Fisher Scientific, Rockford, IL, USA) coupled with a high-performance liquid chromatography system (EASY-nLC 1200, Thermo Fisher). MS raw files generated by LC-MS/MS were searched against the NCBI human Refseq protein database (released on 04-07-2013; 32,015 entries) using MaxQuant (version 1.6.2.10) software enabled with Andromeda search engine. Protease was Trypsin/P. Up to 2 missed cleavages were allowed. Carbamidomethyl (C) was considered as a fixed modification. For the proteome profiling data, variable modifications were oxidation (M) and acetylation (Protein N-term). In this experiment, we performed three repeats of each CRC cell and produced a total of 168 proteome profiles ([330]Figure S4G). During the MS detection, we made a uniformed quality control, and finally resulted in the identification of 6,788 and 7,765 GPs, respectively, in pre-irradiation and post-irradiation CRC cells at a 1% global protein FDR ([331]Table S6). We screened the differentially expressed proteins in pre-irradiation and post-irradiation HSF1-overexpressing (HSF1-OE) or HSF1-knockdown (HSF1-KD) CRC cells (FC > 1.5 or < 0.67). Pathway enrichment analysis was performed according to CPDB database. Western blotting For western blotting, cells were lysed with 0.5% NP-40 buffer containing 20 mM Tris-HCl (pH 8.0), 100 mM NaCl, 1 mM EDTA, 1 mM PMSF, and NONIDET P-40 SUBSTITUTE. The protein concentration was quantified using Bradford assay. For each sample, 30 μg of protein extract was separated using 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis and transferred to nitrocellulose membranes. After blocking with 5% milk (BD Science) solution in TBST (Tris-buffered saline with Tween) for 1–2 h, the membranes were incubated with TBST containing the appropriate primary antibodies overnight at 4°C, followed by a 2 h incubation with horseradish peroxidase-conjugated secondary antibodies. The target protein bands were detected using the Chemiluminescent detection reagent. The mouse monoclonal anti-β-Actin antibody (1:2000, abmart, catalog No: P30002M), and the rabbit Polyclonal anti-HSF1 antibody (1:1000, Signalway Antibody, catalog No: 41037-1) were used, and their specificity was confirmed by western blotting. Western blot quantification was performed using ImageJ software ([332]Version 1.52a, National Institutes of Health, MD, USA). Cytotoxicity assays Cytotoxicity or cell sensitivity to X-rays was analyzed using Cell Counting Kit-8 (CCK8; Beyotime, C0046) assays. For the CCK8 assay, cells were seeded into 96-well plates at a density of 1,000 cells/well. After 24 h, the cells were treated with X-rays using RadSource/RS2000pro-225 at the indicated doses and then cultured for 10 days. CCK8, at a final concentration of 10%, was added to each well and incubated for 1.5 h. Absorbance was measured at 450 nm and 630 nm and the data were analyzed using GraphPad Prism 8 software. Immunofluorescence Cells were irradiated with 2 Gy, collected at the indicated times. Cells in confocal dishes were fixed in 4% paraformaldehyde, permeabilized using 0.5% Triton X-100, and blocked with 3% BSA (Beyotime). The cells were then incubated with anti-γ-H2AX antibody (1:200; ABclonal; AP0687) at 4°C overnight, followed by incubation with FITC-labeled Goat Anti-Rabbit IgG (1:500; Beyotime; A0562) for 1 h at room temperature in the dark. The samples were co-stained with DAPI (Sigma) and examined by confocal laser scanning microscopy. Immunohistochemistry staining and evaluation A standard immunohistochemistry (IHC) protocol was followed to stain the tumor tissue samples using the rabbit polyclonal antibody against, the rabbit polyclonal antibody against HSF1 (Signalway Antibody, catalog No: 41037-1), and the rabbit polyclonal antibody against HDAC6 (Signalway Antibody, catalog No: 40999-1). IHC evaluation was analyzed using an IHC profiler compatible plugin with integrated options for the quantitative analysis of digital IHC images stained for cytoplasmic or nuclear proteins.[333]^65 Moreover, the intensity of the cytoplasmic staining and the percentage of positively stained tumor cells were also scored numerically. Quantification and statistical analysis The proteomic subtypes generated by consensus clustering analysis and their biological function The protein expression matrix of the CRC cohort was used to identify the proteomic subtypes using the consensus clustering method implemented in the R package ConsensusClusterPlus v.3.8.[334]^52 Prior to the consensus clustering analysis, we performed a log[10] transformation and median centered normalization to facilitate the interpretation of the expression data. Then, the top 3,000 proteins with the highest median absolute deviation were subjected to ConsensusClusterPlus in R v.3.5.1 for unsupervised consensus clustering. The cluster analysis was performed with the following setting: maxK = 10, reps = 1,000, pItem = 0.8, pFeature = 1, clusterAlg = "hc", distance = "pearson" for the clustering runs. A preferred cluster result was selected by considering the profiles of the consensus cumulative distribution function (CDF) and delta area under the CDF curve for clustering solutions between 2 and 10 clusters. As shown in [335]Figure S2A, the rank survey profiles of the consensus CDF and the delta area under the CDF curve, along with the consensus membership heatmaps, indicated a four-subtype solution for Fudan CRC cohort using the proteomic data. This showed clear separation and the significant therapy response differences among 4 clusters ([336]Figures S2A and S2B). To generate the abundance heatmap shown in [337]Figure 2A, the CRC samples in each subtype were rearranged from G-I to G-IV, using the signature protein abundance matrix enriched in the signature pathways for each subtype. The signature proteins of each subtype defined here should meet the following criteria: 1) detected in at least 10% of the subtype, 2) differentially expressed of the subtype compared with other subtypes with fold change (FC) > 2 and p < 0.05 (two-sided Wilcoxon rank-sum test). ConsensusPathDB (CPDB) molecular interaction data obtained from 31 different public repositories,[338]^66 and determined the dominant bioprocesses of each subtype. Pathways with a p-value less than 0.05 were regarded to be significant enrichment (two-sided Fisher’s exact test). Correlation between proteomic subtype and clinical characteristics To explore the association of clinical baseline characteristics (including age, gender, response, RECIST, primary site, metastasis, and RAS mutation, etc.) with proteomic subtypes, Fisher’s exact test was performed on categorical variables and one-way ANOVA analysis was applied for comparison for the continuous variables (GraphPad Prism 8 software). p-values less than 0.05 were considered as significantly different. Validation of the fudan subtyping in an independent cohort We validated the Fudan subtyping identified in our cohort using an independent CRC cohort (TCGA cohort) using Rapidminer 9.6.0 (RapidMiner Inc, Boston, USA). Polynomial by Binomial Classification operator uses a binomial classifier and generates binomial classification models for different classes and then aggregates the responses of these binomial classification models for classification of polynomial label. The Fast large margin operator using logistic regression, applied in the subprocess of the Polynomial by Binomial Classification operator, was employed to build the prediction model based on signature proteins of each subtype in the Fudan cohort. The assessment of drug sensitivity The gene expression profiles of CRC cell lines were from Expression 21Q2 Public dataset. These RNASeq files were aligned with STAR and quantified with RSEM, then TPM-normalized and log[2](TPM+1) transformation. Reported values of genes were log[2] (TPM+1). The drug sensitivity (half maximal inhibitory concentration [IC50]) to oxaliplatin (BRD: BRD-K78960041-001-03-2), and 5-fluorouracil (BRD: BRD-K24844714-001-24-5) using Drug sensitivity (PRISM Repurposing Primary Screen) 19Q4 dataset filtered by CRC cell lines. Reported values of drug sensitivity were log[2] [fold change (treatment group vs. control group)]. The Pearson’s correlation of gene expression value [log[2] (TPM+1)] and drug sensitivity [log[2] fold change] were calculated. All datasets were downloaded from DepMap database ([339]https://depmap.org/portal/). Differential protein and pathway analysis in subcohorts The protein expression matrices of the subcohorts were used to perform the differential expression analysis of the sensitive group (S) and the non-sensitive group (NS). The overrepresented proteins of S and NS of subcohorts were defined as proteins differentially expressed in the S and NS groups (NS vs. S > 1.5 or < 0.67); the significantly differentially expressed proteins was defined as proteins with more than 1.5-fold change and Wilcoxon rank-sum test with a Benjamini-Hochberg (BH) adjusted p value cutoff (BH p value < 0.05) ([340]Tables S4, [341]S5, and [342]S7). Pathway enrichment analysis of the overrepresented proteins was performed by CPDB databases. Pathways with a p-value less than 0.05 were regarded to be significant enrichment (two-sided Fisher’s exact test). Gene set enrichment analysis (GSEA) was also used for pathway enrichment analysis.[343]^67 GSEA analysis used by the clusterProfiler R package.[344]^68 GSEA evaluates and determines whether a priori defined sets of genes show statistically significant, cumulative changes in gene expression that are correlated with a specific phenotype, which was applied for pathway analysis in many researches.[345]^69^,[346]^70^,[347]^71 Samples grouped according to S and NS were subjected to GSEA, respectively. Molecular Signatures Database (MSigDB) of hallmark gene sets (H), curated gene sets (C2) and GO gene sets (C5) were used for enrichment analysis. Nominal p values were calculated as Phenotype-based permutation test. A p value of 0.05 was used as a cutoff. The normalized enrichment score (NES) in GSEA was calculated by first ranking the proteins from the most to least significant with respect to S and NS, the entire ranked list was then used to assess how the proteins of each gene set were distributed across the ranked list. Protein-protein interactions The Protein-protein network view summarizes the network of predicted associations for a particular group of proteins by querying the STRING database.[348]^71^,[349]^72 The network nodes are proteins. The edges represent the predicted functional associations. The edges are draw according to the view settings. In evidence mode, an edge may be drawn with up to 7 differently colored lines - these lines represent the existence of the seven types of evidence used in predicting the associations. Red line - indicates the presence of fusion evidence, Green line - neighborhood evidence, Blue line - cooccurrence evidence, Purple line - experimental evidence, Yellow line - text mining evidence, Light blue line - database evidence, Black line - coexpression evidence. In confidence mode the thickness of the line indicates the degree of confidence prediction of the interaction. In addition, the hub proteins were identified through a web source database ([350]https://www.networkanalyst.ca/). The candidate signature proteins were uploaded to the online tool to obtain protein–protein interaction (PPI) information, and the results of the PPI analysis were visualized ([351]Figure 6B). We identified 6 hub proteins (including BAG2, SMYD3, PRKG1, MPHOSPH6, PPFIA1, and NCKIPSD) from 64 signature proteins with criteria of filtering degree > 15 in the S group; and 9 hub proteins (including HDAC6, SMARCA4, RELB, RIPK2, ARL6IP1, SF3B2, CHAF1B, PPP4C, and PPM1B) from 47 signature proteins with criteria of filtering degree > 20 in the NS group. The core hub protein was designated as one had the most interactions with other proteins. The single sample gene set enrichment analysis (ssGSEA) All scores were inferred by single sample gene set enrichment analysis (ssGSEA) method from the GSVA R package based on the protein expression matrix. The genset (c2.all.v7.4.symbols) of Molecular Signature Databse(MSigDB) was used to ssGSEA.[352]^73 The parameters: min.sz = 10, max.sz = 300 were set and other parameters were used default. Construction and validation of predictive models for therapy response Multiple logistic regression analysis was used to construct the therapeutic response prediction model based on the significantly differentially expressed proteins in S and NS groups of subcohorts using in the R software v3.5.1. The backward stepwise method was utilized to feature selection. Samples was randomly divided into 80% of individuals (the training set) and the remaining 20% (the testing set).[353]^74 The 10-flod cross validation and 3-fold cross validation measurements were used to consider probability overfitting on training set, and accuracy, sensitivity and specificity were used to determine the performance of model. For 10-fold cross validation, the data are randomly split into 10 roughly equal-sized parts. One of these parts was held out for validation, and the model was fit on the remaining 9 parts. This fitted model was used to compute the AUC, sensitivity and specificity on the omitted part, and this process was repeated for each of 10 parts. In each iteration among 10 repeats, there were 27, 22 and 11 samples for training model and the remaining 3, 2 and 1–2 samples for internal validation for FOLFOX, FOLFOX plus radiotherapy, and FOLFOX plus cetuximab therapy, respectively. For 3-fold cross validation, the data are randomly split into 3 roughly equal-sized parts. One of these parts was held out for validation, and the model was fit on the remaining 2 parts. This fitted model was used to compute the AUC, sensitivity and specificity on the omitted part, and this process was repeated for each of 3 parts. In each iteration among 3 repeats, there were 20, 16 and 8 samples for training model and the remaining 10, 8 and 4 samples for internal validation for FOLFOX, FOLFOX plus radiotherapy and FOLFOX plus cetuximab therapy, respectively. This process was also conducted using the function of “trainControl” for cross validation in the “caret” package (version 6.0–86) in the R software v3.5.1, of which the parameter of number for cross validation was set as 10 and 3, respectively. The same strategy was adopted for the construction of predictive models in our previously published paper.[354]^60 Moreover, the model diagnostic value was validated in other independent validation cohorts. The predictive models for predicting response of FOLFOX and FOLFOX + cetuximab therapies was separately applied and validated in other public datasets including [355]GSE19860, [356]GSE5851 and Project ID [357]OEP000729 ([358]Figure S5E).[359]^13^,[360]^75 Statistical analysis Standard statistical tests were used to analyze the clinical data, including but not limited to Student’s t test, Wilcoxon rank-sum test, Fisher’s exact test, Pearson’s correlation test, Spearman correlation test, one-way ANOVA. The Wilcoxon rank-sum test was used to examine whether proteins were differentially expressed between the subtype (G-I, G-II, G-III, or G-IV) and other subtypes, S and NS groups, or RAS wildtype and RAS mutation. For correlation analysis, Spearman or Pearson correlation was used. Unless otherwise specified, all statistical tests were two-sided, and statistical significance was considered when p value < 0.05. To account for multiple-testing, the p values were adjusted using the Benjamini-Hochberg FDR correction. All the analyses of clinical data were performed in R (version 3.5.1) and GraphPad Prism 8 software. The p values less than 0.05, 0.01, 0.001, 0.0001 were marked with ∗, ∗∗, ∗∗∗, ∗∗∗∗, respectively. For functional experiments, each was repeated at least three times independently, and results were expressed as mean ± SD. Statistical analysis was performed using GraphPad Prism. This study did not generate custom computer code. No special code was used in this study. Acknowledgments