Prognostic implications of an autophagy-based signature in colorectal cancer.

BACKGROUND: The heterogeneity of colorectal cancer (CRC) poses a significant challenge to the precise treatment of patients. CRC has been divided into 4 consensus molecular subtypes (CMSs) with distinct biological and clinical characteristics, of which CMS4 has the mesenchymal identity and the highest relapse rate. Autophagy plays a vital role in CRC development and therapeutic response.
METHODS: The gene expression profiles collected from 6 datasets were applied to this study. Network analysis was applied to integrate the subtype-specific molecular modalities and autophagy signature to establish an autophagy-based prognostic signature for CRC (APSCRC).
RESULTS: Network analysis revealed that 6 prognostic autophagy genes (VAMP7, DLC1, FKBP1B, PEA15, PEX14, and DNAJB1) predominantly regulated the mesenchymal modalities of CRC. The APSCRC was constructed by these 6 core genes and applied for risk calculation. Patients were divided into high- and low-risk groups based on APSCRC score in all cohorts. Patients within the high-risk group showed an unfavorable prognosis. In multivariate analysis, the APSCRC remained an independent predictor of prognosis. Moreover, the APSCRC achieved higher prognostic power than commercialized multigene signatures.
CONCLUSIONS: We proposed and validated an autophagy-based signature, which is a promising prognostic biomarker of CRC patients. Further prospective studies are warranted to test and validate its efficiency for clinical application.

Introduction

Colorectal cancer (CRC) is one of the most common malignancies in the world and a leading cause of cancer-related death.[ In 2017, the prevalence rate of CRC in China ranked the third among all malignant tumors, about 37.6 per 100,000, and the mortality rate ranked the fourth, about 19.1 per 100,000.[ Although great progress has been made in CRC treatment including surgery and adjuvant chemotherapy, the 5-year survival of CRC patients is still not optimistic due to the high recurrence rate.[ At present, tumor staging remains the gold standard for estimating the prognosis of CRC. However, the clinical prognosis of many CRC patients with the same stage or treatment strategy varies.[ The genetic heterogeneity of patients greatly affects intrinsic clinical diversity.[ Therefore, there is a need to integrate tumor heterogeneity to identify novel prognostic predictors for CRC.

Over the past 20 years, researchers have identified a number of CRC prognostic marker genes and built multigene-based prognostic signatures that can assess the prognostic risk of CRC.[ Due to the inherent genetic heterogeneity of CRC, the robustness of most markers is less than expected. In 2015, CRC was classified into 4 consensus molecular subtypes (CMS) with distinct molecular and clinical characteristics, among which CMS4 has mesenchymal characteristics.[ CMS2 subtype responded well to anti-EGFR and HER2 inhibitors, while CMS1 and CMS4 showed higher response rates to HSP-like inhibitors. The CMS4 subtype showed significant chemoresistance to the combination of 5-Fu and Luminespib.[ Importantly, the CMS4 shows the highest recurrence rate and the worst relapse-free survival (RFS) indicating the need to integrate the intrinsic modalities of this malignant subtype for risk management in CMS4.

Autophagy is an important mechanism to degrade cytoplasmic components and maintain intracellular stability[ in many diseases, especially cancer.[ Autophagy is thought to suppress early-stage tumor development and promote late-stage tumor development.[ Many studies have explored the mechanism of autophagy-related genes (ATGs) in the development and drug-resistance of CRC.[ Therefore, ATGs is expected to be a new prognostic indicator for CRC.

Considering the highly heterogeneous nature of CRC, we integrate molecular modalities and ATGs underlying CMS4. We pooled and analyzed 6 public cohorts containing 1428 CRC patients to develop and validate an autophagy gene-based prognostic signature for CRC (APSCRC). Although autophagy prognostic markers for CRC have already been mentioned before,[ risk assessment incorporating molecular subtypes has not been investigated. Our signature incorporates the autophagy system and tumor heterogeneity and would be useful for selecting patients with CRC to give more precise treatment.

Materials and methods

Ethics statement

The researchers were authorized to conduct the study by the Ethics Committee of the Beilun People's Hospital, Ningbo, China. All procedures were implemented in accordance with the Declaration of Helsinki and relevant policies in China.

Patient datasets

In this study, we collected gene expression profiles (GEPs) from 6 independent datasets (Supplemental Table 1), including 1428 cases. These datasets involved patients from the GSE14333 (n = 290),[ GSE17538 (n = 232),[ the GSE39582 (n = 564),[ the GSE33113 (n = 90),[ the GSE38832 (n = 122),[ and the GSE37892 (n = 130).[ Expression data and corresponding clinical information for all datasets were downloaded from Gene Expression Omnibus. Molecular subtyping information was retrieved from Guinney's study.[ Detailed clinical characteristics of these 6 datasets were described in Table 1. The design and workflow of this study were shown in Figure 1A. For each cohort, the GEPs probe IDs were transformed to genes symbols, if multiple probe IDs correspond to the same genes symbol, the one with the highest mean value was kept as the representative of the corresponding gene.[

Table 1

Patients’ characteristics in public data sets.

	Training cohort	Validation cohorts
	GSE14333 (n = 290)	GSE17538 (n = 232)	GSE39582 (n = 564)	GSE33113 (n = 90)	GSE37892 (n = 130)	GSE38832 (n = 122)
Age (y)	66 (26–92)	65 (23–94)	67 (22–97)	70 (34–95)	68 (22–97)
Gender
Female	126 (43%)	110 (47%)	256 (45%)	48 (53%)	61 (47%)
Male	164 (57%)	122 (53%)	309 (55%)	42 (47%)	69 (53%)
Location
Left	161 (56%)		342 (61%)		72 (55%)
Right	128 (44%)		223 (39%)		57 (44%)
Unknown					1 (1%)
Stage
I	44 (15%)	28 (12%)	32 (6%)			18 (15%)
II	94 (32%)	72 (31%)	264 (47%)	90 (100%)	73 (56%)	35 (29%)
III	91 (31%)	76 (33%)	205 (36%)		57 (44%)	39 (32%)
IV	61 (21%)	56 (24%)	60 (11%)			30 (25%)
Unknown			4 (1%)
MSI status
msi			75 (13%)	25 (28%)
mss			443 (78%)	65 (72%)
Unknown			47 (8%)
BRAF
Wild			460 (81%)	73 (81%)
Mut			51 (9%)	17 (19%)
Unknown			54 (10%)
Subtype
CMS1	63 (22%)	43 (19%)	91 (16%)	20 (22%)	11 (8%)	20 (16%)
CMS2	119 (41%)	72 (31%)	232 (41%)	34 (38%)	49 (38%)	29 (24%)
CMS3	50 (17%)	36 (16%)	69 (12%)	10 (11%)	19 (15%)	19 (16%)
CMS4	58 (20%)	47 (20%)	126 (22%)	21 (23%)	39 (30%)	33 (27%)
Unknown		34 (15%)	47 (8%)	5 (6%)	12 (9%)	21 (17%)

CMSs = consensus molecular subtypes.

Figure 1

Network inference identifies 6 autophagy genes as key regulators of the CMS4 subtype of CRC. Study design of the current work (A). Integrated network displays the relationships between the 6 autophagy genes (VAMP7, DLC1, FKBP1B, PEA15, PEX14, and DNAJB1) and target genes expression data (B). CMSs = consensus molecular subtypes, CRC = colorectal cancer.

Integrated network analysis

ATGs were downloaded from the HADb database.[ ATGs measured by all datasets were screened. Regulatory network inference was applied to investigate the relationship between ATGs and potential target genes underlying CMS4. Together, we used the GSE14333 dataset as the training cohort. Forty-six ATGs (|log2FC| > 0.25, P < .05) and 1407 target genes (|log2FC| > 0.5, P < .05) were determined differentially expressed in CMS4 compared with the other 3 subtypes (CMS1, CMS2, and CMS3). Integrated analysis was performed by the “RTN” (R package, version 2.10.0).[ Master regulator analysis (MRA) was performed to examine the overrepresentation of epithelial–mesenchymal transition (EMT) signature[ within each autophagy gene's regulon.

Development of the autophagy-based prognostic signature for CRC (APSCRC)

Six autophagy genes were the master regulatory factors and differentially expressed in CMS4 which is the worst survival molecular subtype in CRC. Based on these 6 autophagy genes, the multivariable cox proportional-hazards model was used to construct the APSCRC in the training cohort. The risk-score formula was constructed as follows: risk score = (0.1991 x DLC1) + (0.2351 x FKBP1B) + (0.0973 x PEA15) + (0.5989 x DNAJB1) +(−0.1977 x VAMP7) +(−0.1584 x PEX14).

Validation of the APSCRC

The prognostic relevance of APSCRC was further evaluated in 5 independent validation cohorts by Kaplan–Meier analysis, respectively. To test whether the APSCRC can be considered as an independent prognostic predictor of CRC patients, univariate and multivariate cox analyses were performed with other clinical factors. In the multivariable cox regression, gender, tumor locations, microsatellite instability (MSI) status, and BRAF mutation status were included as covariates.

Gene set enrichment analysis (GSEA)

GSEA[ was conducted using the “fgsea” package (Bioconductor package, version 1.12.0) with 1000 permutations. Gene sets were retrieved from the Molecular Signature Database (MSigDB hallmark and kegg, version 7).[ A P value below .05 was used to choose significant gene sets.

Statistical analysis

Differential analysis of autophagy genes and target genes between CMS4 subtype and non-CMS4 subtypes was done by package “limma.”[ For each cohort, the upper quantile risk value was used to stratify patients into high- and low-risk groups. Kaplan–Meier analysis was performed using the log-rank test using “survival” (R package, version 2.41.3). Univariate cox regression analysis was conducted to show the prognostic value for the selected autophagy signature. Uni- and multivariable analyses were conducted by the Cox proportional hazards model. The subtype prediction power of APSCRC was determined by receiver operating characteristic curves. For all tests, a P value below.05 was used to choose significant gene sets. Statistical significance is presented as following ∗P < .05, ∗∗P < .01, ∗∗∗P < .001. All the statistical tests were conducted using R (version 3.6.1).

Results

Integrative analysis reveals 6 autophagy genes as key regulators in the CMS4 subtype of CRC

CRC is a heterogeneous disease. In recent studies,[ 4 CMSs have been identified, among which CMS4 has the worst RFS and overall survival (OS) (Supplemental Fig. 1A, B). Therefore, we intended to integrate the molecular modalities under this subtype to improve risk assessment of CRC hereafter. Focusing on the CMS4 subtype of CRC, we applied a network-based approach to investigate the regulatory role of the autophagy system underlying this particular subtype. A total of 1428 CRC patients from 6 independent public cohorts were included (Table 1). Two hundred thirty-two autophagy-related genes (ATGs) were downloaded from the HADb database.[ Based on the GSE14333 dataset, we performed differential analysis. As a result, 46 ATGs (|log2FC| > 0.25, P < .05) and 1407 target genes (|log2FC| > 0.5, P < .05) were differentially expressed in the CMS4 subtype compared with non-CMS4 subtypes (Fig. 1A). Under the expression profiles of these prioritized autophagy genes and target genes, a regulatory network was constructed by calculating the mutual information between an autophagy signature and its potential targets[ (Fig. 1B). MRA identified 6 autophagy genes (VAMP7, DLC1, FKBP1B, PEA15, PEX14, and DNAJB1) as the key factors of CMS4 (Supplemental Table 2). Univariable Cox regression analysis demonstrated prognostic association of the 6 autophagy genes with RFS in more than 1 set of cohorts (Fig. 2A–D). Therefore, the network-based approach identified 6 subtype-specific autophagy genes, with prognostic value in CRC, as key regulators of the CMS4 subtype of CRC and can be potentially applied for risk assessment of CRC patients.

Figure 2

Forest plots show the prognosis association of these 6 autophagy genes in GSE14333 (A), GSE39582 (B), GSE17538 (C), and GSE33113 (D). (∗P < .05, ∗∗P < .01, ∗∗∗P < .001).

Development of an autophagy-based prognostic signature for CRC (APSCRC)

Using the GSE14333 cohort as the training set, a risk model designated APSCRC was constructed with these 6 autophagy genes as following formula: risk score = (0.1991 x DLC1) + (0.2351 x FKBP1B) + (0.0973 x PEA15) + (0.5989 x DNAJB1) +(−0.1977 x VAMP7) +(−0.1584 x PEX14). Subsequently, risk scores were calculated for all the training and validation cohorts (Supplemental Table 3). The risk scores showed high efficiency in detecting the CMS4 subtype, suggesting that the APSCRC correlated with mesenchymal phenotype (Fig. 3A–D). The upper quartile risk value was set as the cutoff to separate patients into high- and low-risk groups across all datasets. In the training set, the high-risk group patients displayed significantly poorer RFS than patients in the low-risk group (HR: 2.71, 95% CI: 1.55 - 4.77; P = 2.91 x 10-4) (Fig. 4A). When considering patients with stage II CRC only, the APSCRC remained prognostic in terms of RFS for the training set (HR: 3.77, 95% CI: 1.31-10.9; P = 8.44 x 10-3) (Fig. 5A).

Figure 3

ROC curves show the detection efficiency of the CMS4 subtype for APSCRC in GSE39582 (A), GSE17538 (B), GSE33113 (C), and GSE37892 (D). CMSs = consensus molecular subtypes, ROC = receiver operating characteristic.

Figure 4

Prognostic value of the APSCRC in CRC. Kaplan–Meier survival analysis showed that the high-risk group has an unfavorable RFS in GSE14333 (A), GSE17538 (B), GSE39582 (C), GSE33113 (D), GSE37892 (E), GSE38832 (F), and the pooled validation cohorts (G). P values are calculated by the log-rank tests. APSCRC = autophagy-based prognostic signature for CRC, CRC = colorectal cancer.

Figure 5

Prognostic value of the APSCRC for stage II CRC within GSE14333 (A), GSE17538 (B), GSE33113 (C), and the pooled validation cohorts (D). P values are calculated by the log-rank tests. APSCRC = autophagy-based prognostic signature for CRC, CRC = colorectal cancer.

Validation of the APSCRC for survival prediction

To verify the prognostic power of the APSCRC, we calculated the survival difference within the 2 risk groups in 5 validation cohorts. As expected, the high-risk group showed significantly reduced RFS (P =.03 to P = 5.04 x 10−5) (Fig. 4, B–F) and OS (P =.003 to P = 7.62 x 10−6) (Fig. 6A–C) compared with the low-risk group in 5validation cohorts. When considering patients within stage II CRC only, the APSCRC remained prognostic in terms of RFS for validation cohorts (P = .04 to P = .006) (Fig. 5, B and C). When pooled all validation datasets together, a more obvious RFS (HR: 1.96, 95% CI: 1.58-2.42; P < 3.07 x 10−10) (Fig. 4G), stage II RFS (HR: 2.12, 95% CI: 1.44-3.12; P = 1.02 x 10−4) (Fig. 5D) and OS (HR: 1.94, 95% CI: 1.54-2.45; P = 1.52 x 10−8) (Fig. 6D) differences were observed between the high- and low-risk groups. Moreover, the APSCRC also showed prognostic value for stage I, III, and IV CRC patients, in addition to stage II CRC (Supplemental Fig. 2A–C). To testify whether the APSCRC was an independent prognostic predictor, uni- and multivariate Cox regression analyses were conducted. After adjusting for gender, tumor locations, MSI status, and BRAF mutation status, the APSCRC remains an independent predictor of prognosis (Table 2).

Figure 6

Kaplan–Meier plots showing differences in overall survival among different risk groups within GSE17538 (A), GSE39582 (B), GSE38832 (C), and the pooled validation cohorts (D). P values are calculated by the log-rank tests.

Table 2

Univariate and multivariate analysis of autophagy signature and clinicopathological factors.

	Pooled validation cohorts
	Univariate		Multivariate
	HR (95% CI)	P	HR (95% CI)	P
Gender (male vs female)	1.12 (0.89–1.42)	.33	1.07 (0.76–1.51)	.68
Location (left vs right)	1.22 (0.92–1.61)	.16	1.36 (0.94–1.97)	.11
MSI (msi vs mss)	2.48 (1.44–4.28)	1.00E-03	4.54 (2.08–9.90)	1.0E-04
BRAF (mut and wild)	1.01 (0.62–1.65)	.97	1.98 (0.99–3.99)	.05
APSCRC (high vs low risk)	1.79 (1.40–2.86)	3.50E-06	2.00 (1.35–2.97)	5.0E-04

APSCRC = autophagy-based prognostic signature for CRC.

Biological annotation of the APSCRC

To gain insight into the biological differences between risk groups, we performed GSEA analysis. Enrichment analysis between high- and low-risk groups identified that many mesenchymal phenotype-related pathways, including the TGF-beta signaling, epithelial–mesenchymal transition, and focal adhesion, were positively enriched in the high-risk group (Fig. 7A–D and Supplemental Table 4). When compared with a clinically applicable and commercialized biomarker Oncotype Dx Colon Cancer,[ the APSCRC achieved a higher C-index in training and validation cohorts (Fig. 8A and B).

Figure 7

The enrichment of dysregulated pathways between the high- and low-risk groups.

Figure 8

The comparison of C-index between APSCRC and Oncotype DX Colon Cancer (A). C-index values in different cohorts (B). APSCRC = autophagy-based prognostic signature for CRC.

Discussion

CRC is one of the most common types of cancer worldwide and the fourth leading cause of cancer-related deaths.[ In recent years, researchers have constructed numerous multigene-based prognostic biomarkers that can divide CRC patients into different risk groups.[ However, due to the genetic heterogeneity of CRC, the robustness of most markers is less than expected. A deeper understanding of the molecular characteristics of CRC may lead to better risk stratification of CRC patients. In 2015, CRC was classified into 4 CMS with distinct molecular and clinical characteristics.[ Patients within CMS4 had the highest recurrence rate and the worst RFS. Autophagy is thought to suppress early-stage tumor development and promote late-stage tumor development.[ Autophagy can be classified into 3 major types; macro-autophagy, micro-autophagy, and chaperone-mediated autophagy.[ Inhibition of macro-autophagy promotes CMA-dependent degradation of mutant p53.[ Higher expression of LAMP-2A indicated activated CMA in CRC.[ However, these processes have not been well studied in CRC. Micro-autophagy is a direct engulfment of cellular components through late indentation of the late endosome. There is evidence in lung cancer that amino acid starvation, an important factor in cancer growth, causes microautophagy.[ However, the role of micro-autophagy in CRC has not been studied in detail. Many studies have explored the mechanism of ATGs in the development and drug resistance of CRC.[ Meanwhile, autophagy-based prognostic markers for CRC have also been reported.[ However, most studies only focus on autophagy without considering the heterogeneity of CRC. To our knowledge, no research has been done for risk stratification by integrating the autophagy system and the mesenchymal modalities of CMS4 in CRC.

In this study, we identified 6 autophagy genes as the key factors of CMS4 and constructed an APSCRC. Within these 6 autophagy genes, FKBP1B is involved in an autophagy-based prognostic panel for multiple myeloma.[ PEA15 shows a critical role in chemoresistance and metastasis in CRC. Meanwhile, the PEA15 expression level is significantly associated with the prognosis of CRC patients.[ DNAJB1 actively regulates the EGFR signaling pathway by destabilizing MIG6.[ DLC1 is a direct target gene for miR-106b and correlates with prolonged survival in CRC.[ Each patient was assigned a risk score by the APSCRC in all cohorts. Patients within each cohort could be separated into high- and low-risk groups based on the upper quartile risk value. The APSCRC was proved as a stable prognostic biomarker of CRC, and patients with high risk scores were correlated with poor RFS and OS in 5 independent validation cohorts. Cox regression analyses indicated that the APSCRC was an independent prognostic predictor after adjusting several clinical risk factors. Moreover, some mesenchymal phenotype-related pathways, such as EMT, TGF-beta, and focal adhesion, were positively enriched in the high-risk group. When compared with a clinically applicable and commercialized biomarker Oncotype Dx Colon Cancer,[ the APSCRC achieved a higher C-index. Our autophagy gene-based signature can effectively predict patient survival of CRC patients.

This study still has some limitations. First, the prognostic signature was screened from gene expression profiles generated from microarray platforms, which are expensive, difficult to operate, and involve professional bioinformatics expertise, so it is difficult to be popularized in daily clinical application. Second, the training and validation data sets were all from retrospective studies in the study. In the following improvement process, more datasets containing more clinical characteristics need to be included for more extensive screening and validation. Finally, a uniform risk threshold needs to be established to facilitate potential clinical applications.

Taken together, we applied a network-based approach to construct an autophagy-based prognostic signature for CRC. Our study is the first attempt to integrate colorectal cancer heterogeneity and the autophagy system to develop the prognostic signature for CRC.

Author contributions

Conceptualization: peng shu, Xingguo Zhang.

Data curation: Liangbin Wang, Xinlei Jiang.

Formal analysis: Liangbin Wang, Xinlei Jiang.

Investigation: Xingguo Zhang.

Methodology: Liangbin Wang.

Project administration: peng shu, Xingguo Zhang.

Software: Liangbin Wang.

Supervision: peng shu.

Validation: Liangbin Wang.

Visualization: Liangbin Wang.

Writing – original draft: peng shu, Liangbin Wang.

Writing – review & editing: peng shu, Xingguo Zhang.