The transcriptome profile, mutation data and clinical follow-up information of 414 samples were downloaded from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). A prospectively collected colon cancer cohort from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) was consistent of paired tumor and non-tumor colon tissues from 110 colon cancer patients (18). The RNA sequencing data, mutation data and clinical data were obtained from the cBioPortal database (https://www.cbioportal.org/study/summary?id=coad_cptac_2019). The two RNA-sequencing cohorts were merged into one combined RNA-seq cohort after removing the batch effect. A total of 261 TP53-mutant samples in the combined RNA-seq cohort were used for further analysis. The data we used was obtained from public database and the patients involved in the database have obtained ethical approval.

A total of three microarray datasets were enrolled in this study. GSE39582 and GSE17536 were retrieved from the GPL570 platform, and GSE41258 were retrieved from the GPL96 platform. The expression data and corresponding clinical information of 585 samples in the GSE39582, 177 samples in the GSE17536 and 186 samples in the GSE41258 were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). After eliminating samples without complete clinical data and survival data, 541 samples in the GSE39582, 177 samples in the GSE17536, and 182 samples in the GSE41258 were used for further analysis. Notably, 189 samples in the GSE39582 were recorded as TP53 mutate while other 352 samples were recorded as TP53 wild-type, and the TP53-mutant samples were used for signature establishment. GSE17536 and GSE41258 were combined into one GEO cohort as validation.

One hundred and eighty-nine samples with TP53 mutations in the GSE39582 were randomized into the training set (N = 94) and the testing set (N = 95). LASSO-Cox regression method (19, 20) was performed by the “glmnet” R package in the training set and the risk score of each sample were calculated based on the coefficient and expression level of each selected genes.

To validate our result in a larger sample size resulting in achieving a more reliable finding, we performed RF analysis via the “RandomForest” R package to construct a TP53 mutation status prediction model in GSE39582. The TP53 mutation status of samples in the combined GEO cohort was estimated by the TP53 mutation status prediction model. Finally, one hundred and ninety-six samples in the combined GEO cohort were predicted as TP53 mutant samples.

The Broad Institute Cancer Cell Line Encyclopedia (CCLE) project contains the expression profile data and somatic mutation data of over a thousand human cancer cell lines (21) and the data were downloaded from the CCLE website (https://portals.broadinstitute.org/ccle/). CERES score is a computational method to estimate gene dependency levels from CRISPR-Cas9 essentiality screens while accounting for the copy-number-specific effect (22). A more negative CERES score indicates that the gene is essential for cell viability in the certain cell line. The CERES scores of genome-scale CRISPR knockout screens for 18,333 genes in 739 cell lines were acquired from the dependency map (DepMap) portal (https://depmap.org/portal/).

Genomics of Drug Sensitivity in Cancer (GDSC) is the largest public database that provides drug response data of over five hundred compounds in about one thousand human cancer cell lines (23). In the GDSC database, forty-nine colon and rectal adenocarcinomas (COAD-READ) cell lines were found and twenty-five of them harbor TP53 mutations. To analyze the drug response in TP53-mutant colon cancer, we downloaded the IC50 data of the drugs in the TP53-mutant COAD-READ cell lines in the GDSC database (https://www.cancerrxgene.org/). The information on the cell lines was downloaded from the DepMap portal (https://depmap.org/portal/).

The survival curves were drawn by using Kaplan–Meier survival analysis via the “Survival” R package. The difference in survival time between the two groups was tested by a log-rank test. The time-dependent receiver operating characteristic (ROC) curve was performed via the “timeROC” R package, and the value of the area under the curve (AUC) was obtained as the criterion for the accuracy of the prognostic signature. Multivariate Cox analysis was used to evaluate the independence of our prognostic signature on the prognosis in TP53-mutant COAD. Nomogram and calibration curve were fabricated by “rms” R package. All the statistical analysis and visualization were performed by the R (version 4.0.2). A two-tail p-value < 0.05 was considered a significant difference. A flow chart of our study was shown in Figure 1.

The clinical characteristics of all samples enrolled in this study were summarized in the Supplementary Table S1. Firstly, two hundred and fifty-three genes were identified as prognostic genes in TP53-mutant samples in the GSE39582 by univariate cox analysis (Supplementary Table S2, all p < 0.01). And then, TP53-mutant samples in the GSE39582 were randomized into the training set and the testing set. Based on the prognostic genes, we established a 16-gene prognostic signature by LASSO-Cox regression in the training set (Figure 2A). The selected 16 genes included GALK1, TGIF2, TAPBPL, SPINK1, ZNF500, LAMC1, MICB, RPL8, EEF1D, MAPKBP1, ZNF250, RFX3, ETV1, SERINC3, DIP2C, and AKT3. The risk score of each sample was calculated based on the expression level and the coefficient of each gene using the uniform formula: riskscore=∑i=0nexpri∗coefi (Supplementary Table S3). The expression profiles of these sixteen signature genes in GSE39582 were shown in Figure 2B.

The risk scores of each sample were calculated by uniform formula, and the dataset was divided into high and low-risk groups based on the median value of risk scores. The clinical characteristics of each subgroup were summarized in Supplementary Table S4. Samples with high-risk scores had inferior outcomes compared to those with the low-risk score in the training set (Figure 2C, p < 0.0001), the testing set (Figure 2D, p = 0.00077), and the TP53 mutant samples in the GSE39582 cohort (Figure 3A, p < 0.0001). Besides, the AUC values of 1-year, 3-year, and 5-year survival prediction in TP53-mutant samples of GSE39582 were 0.8028, 0.7852, and 0.7557, respectively (Figure 3B). Moreover, the TP53-mutant samples in the combined RNA-seq cohort and the estimated TP53-mutant samples in the combined GEO cohort were enrolled as validation. Similarly, a high-risk score was associated with a worse prognosis in the TP53-mutant samples in the combined RNA-seq cohort (Figure 3C, p = 0.017) and the estimated TP53-mutant samples in the combined GEO cohort (Figure 3D, p = 0.014). In addition, the univariate-cox analysis showed that risk score and M category were prognostic factors in TP53-mutant COAD of GSE39582 (Figure 3E, all p < 0.001), thus we enrolled these two indexes for multivariate-cox analysis. Notably, risk score and M category were both independent poor factors for the prognosis in TP53-mutant COAD (Figure 3E, all p < 0.001). In the combined RNA-seq cohort, risk score, age, T category, N category, and M category were all poor factors for the prognosis in TP53-mutant COAD (Supplementary Figure S1A, all p < 0.05), besides, age and M stage were independent poor factors on the prognosis in TP53-mutant COAD (Supplementary Figure S1A, all p < 0.05). Finally, we fabricated the nomogram based on our signature and multiple clinical features for TP53 mutant patients with microarray data (Figure 3F) or RNA sequencing data (Supplementary Figure S1B). Besides, we also performed the calibration curve to test the predictive efficiency of the nomogram. The result indicated that both two nomograms had good efficiency in the prediction of 1-year to 3-year survival in TP53 mutant COAD patients (Figure 3G, Supplementary Figure S1C).

To test whether the established signature can be used for prognosis prediction in COAD with TP53 wild-type, the prognostic significance of our signature were also evaluated in COAD samples with TP53 wild type, including 352 TP53 wild-type samples in the GSE39582 and 253 samples with TP53 wild-type in the combined RNA-seq cohort. The risk score of each sample were also calculated with uniform formula and each cohort were divided into high risk and low risk groups based on the median value of risk score, respectively. Comparison of survival curve between high risk and low risk groups showed that no significant difference was found in the survival between high risk group and low risk group in both the TP53 wild-type samples in the GSE39582 (Figure 4A, p = 0.12) and samples with TP53 wild-type in the combined RNA-seq cohort (Figure 4B, p = 0.31). Moreover, we estimated the predictive efficiency on the 1-year, 3-year, and 5-year survival of COAD with TP53 wild type in two cohorts using ROC curve analysis. The result both indicated that the established signature had poor predictive efficiency on the prognosis in TP53 wild-type COAD (Figures 4C,D). Taken together, the established signature failed to use for prognosis management in COAD with wild-type TP53, which indicated that our signature was especially used for prognosis prediction in TP53 mutant COAD patients.

The expression data of 49 COAD cell lines were obtained from the CCLE database. Firstly, the correlation between the expression level of genes and risk score was performed, and the expression level of seven genes was positively correlated with a risk score, including LUM, MIA3, NPTN, PDGFRB, PHOQ, SGPP1, and VAMP4 (Figure 5A). And then, we found that 15 of the 49 COAD cell lines harbored TP53 mutations, thus, we further analyzed the relationship between the CERES score of seven above-mentioned genes and the risk score in the 15 TP53-mutant COAD cell lines. The result showed that the risk score was negatively correlated with the CERES score of all seven genes (Figures 5B–H). Finally, we defined three of the seven genes as an essential target for high-risk TP53-mutant COAD, which with CERES scores less than zero in 80% TP53 mutant COAD cell lines, including SGPP1, RHOQ, and PDGFRB.

The sensitivity data of 518 anti-tumor drugs was obtained from the GDSC database which contained two datasets (GDSC1 and GDSC2). We identified candidate drugs by analyzing the IC50 of drugs in both two datasets. We found that the 15 TP53-mutant COAD cell lines were sensitive to eight drugs in the GDSC1 dataset and five drugs in the GDSC2 dataset, including Mitomycin-C, Dacinostat, Belinostat, GSK1059615, Apitolisib, AZD5438, AZD7762, IGFR-3801, Staurosporine, MK-1775, Dinaciclib, Sabutoclax, and MG-132 (Figure 6A). Moreover, we compared the IC50 between samples with high risk and low risk, and the result showed that the high-risk group had lower IC50 compared to the low-risk group in three drugs including IGFR-3801, Staurosporine, and Sabutoclax (Figure 6B), which indicated that IGFR-3801, Staurosporine, and Sabutoclax were potential agents for high-risk COAD with TP53 mutations.