Development and Validation of an Autophagy-Stroma-Based Microenvironment Gene Signature for Risk Stratification in Colorectal Cancer

Background

Colorectal cancer (CRC) is the fourth most common cancer (8.2% of all new cancer cases) and the second leading cause of cancer-related death (8.8% of all cancer deaths) in the USA.¹ However, its incidence varies widely from region to region. CRC can be regarded as a sign of socio-economic development because its incidence rate tends to rise consistently with increase of the human development index.² Therefore, it is foreseeable that its incidence will increase significantly with further development of the global economic situation. It is worth noting that the 5-year relative survival rate for colorectal cancer patients is 65%. A previous study reported that the incidence rate of patients diagnosed with stage I or stage II CRC accounted for 20% and 22%, respectively, and the 5-year relative survival rates were 91% and 82%, respectively. However, the 5-year survival rate for stage IV CRC decreases to 12%.³ In recent years, the tumor microenvironment has become increasingly important to researchers. Studies have shown that the tumor microenvironment of CRC is a key factor affecting the prognosis of colorectal cancer,^4,5 and holds keys to future cancer therapy.⁶ Generally, colorectal tumors are rich in stroma.⁷ Previous studies have reported that stromal cells are differentially expressed in CRC tissues and adjacent tissues,⁸ and affect the occurrence of colorectal cancer.⁹ Notably, stromal cells affect the prognosis of CRC patients by affecting the tumor microenvironment.¹⁰ On the other hand, the role of autophagy in tumors is controversial. Some researchers believe that autophagy can promote the occurrence and development of colorectal cancer.^11–13 However, a recent study by Nassour et al¹⁴ reported that autophagy has an inhibitory effect on cancer. For example, one study has shown that autophagy can inhibit the occurrence and development of colorectal cancer.¹⁵ Interestingly, accumulating evidence indicates that there is a direct or indirect interaction between stromal cells and autophagy in the tumor microenvironment.^16,17 In this study, we speculate that the interaction between autophagy and stromal cells can provide prognostic value for patients with colorectal cancer. Therefore, we combined the status of autophagy and stromal cell into the existing clinicopathological characteristics and staging system through a series of systematic analysis, and developed a new prognostic model aimed at improving the prognosis of colorectal cancer.

Methods

Patient Cohort

The discovery cohort, containing 574 CRC patients, used microarray gene expression data retrieved from the NCBI Gene Expression Omnibus (GSE39582). We then conducted external validations using four independent cohorts, which included 200 patients from the GSE17538 dataset, 92 patients from the GSE38832 dataset, 179 patients from the TCGA database, and 82 colorectal cancer patients from the First Affiliated Hospital of China Medical University (CMU). The 82 CRC patients from the hospital signed informed consent before the study begun and the procedure was approved by the Ethics Committee of The First Affiliated Hospital of China Medical University. Committee approval number: AF-SOP-07-1.1–01

Identification of Autophagy Status and Autophagy-Related DEGs

We used Rtsne¹⁸ R package to randomly simulate the autophagy state. T-SNE utilizes an unsupervised, non-parametric approach to divide patients into several different groups based on given characteristics. We downloaded autophagy marker gene set from the molecular signature database (MSigDB^19,20 version7.1), which included 152 genes. Two groups of autophagy states were determined depending on the clustering results, including the “autophagy^type1” group and the “autophagy^type2” group. In addition, we analyzed the expression changes of target genes involved in autophagy signaling pathway between the “autophagy^type1” group and the “autophagy^type2” group to investigate their association with autophagy. These targets were downloaded from the molecular signature database (MSigDB^19,20 version7.1), including 120 genes associated with “positive regulation of autophagy” and 84 genes associated with “negative regulation of autophagy”. We then used the limma²¹ algorithm to determine differentially expressed genes (DEGs) between two groups.²⁰ Genes with a false discovery rate (FDR)-corrected p-value (q-value) of < 0.0001 and an absolute log2 fold change (log FC) ≥ 1.5 were considered as autophagy-related genes.

Identification of Stroma Status and Stroma-Related DEGs

xCell²² is a tool that performs cell type enrichment analysis from gene expression data for 64 immune and stromal cell types. In this study, we calculated the stromal cell score using “xCell” R package, and divided patients into two groups according to the stromal cell score. Moreover, we examined and compared the prognosis of both groups. Patients with high stromal cell score were classified into the “stroma^high” group based on the median, while those with low stromal cell scores were considered to be in the “stroma^low” group. Next, we identified DEGs between the low and high stroma groups using limma algorithm. Genes with a false discovery rate (FDR)-corrected p-value (q-value) of < 0.0001 and an absolute log2 fold change (log FC) ≥ 1.5 were considered as stroma-related genes.

Identification of Autophagy-Stroma-Related Prognostic DEGs

The above autophagy and stromal states were further composed into a two-dimensional index, and the patients were segmented into three groups, namely autophagy^type1/stroma^low, autophagy^type2/stroma^high, and the mixed group. Autophagy-stroma related DEGs (|log₂ FC|>1.5, FDR correction p<0.0001) were obtained by comparing the expression of autophagy^type1/stroma^low and autophagy^type2/stroma^high group. Next, the protective gene set and risk gene set were identified by overlapping the autophagy-stroma-related genes and autophagy/stroma associated genes. Finally, univariate Cox regression analysis was performed for all protective and risk genes in order to obtain prognostic indicators associated with autophagy and stroma. P <0.05 was considered to be statistically significant.

Derivation of Autophagy-Stroma-Based Gene Classifier and Prognosis Classifier

A popular penalization method is the Least Absolute Shrinkage and Selection Operator (LASSO),^23,24 which can avoid overfitting and retain valuable variables. LASSO has been extended and is used in Cox regression survival analysis and is ideal for high-dimensional data. In our study, we screened the most prognostic gene marker using the LASSO-Cox regression model of glmnet²⁵ R package from identified prognostic factors associated with autophagy and stroma in the discovery cohort. To reduce the potential instability of the results, three cross-validation tests and 1000 iterations were performed. Notably, the best tuning parameter is determined using the 1-SE (standard error) rule. Afterwards, we built a prognostic signature on an individual-level predictive score based on selected risk gene characteristics. The risk score for an individual was the sum of the product of the expression level of the risk marker gene and the coefficient obtained by the LASSO model, that is, the risk score =∑ (coefficient × expression of the characteristic gene). Based on the risk score at the individual-level, an optimal critical value was determined using the maximum selectable rank statistics method²⁶ in order to develop a prognostic classifier for colorectal cancer patients.

Pan-Cancer Data

ENCORI²⁷ was used to analyze gene expression RNAseq dataset of pan-cancer data in TCGA.

Statistical Analysis

Statistical analyses were conducted using R version 3.6.3 (http://www.R-project.org). R package “Rtsne” was used to implement the t-SNE algorithm based on nonlinear dimensionality reduction. Moreover, the “xCell” R package was used to score tumor microenvironment stromal cells, while the “glmnet” R package was used for LASSO-Cox regression modeling. The datum was analyzed using appropriate standard statistical tests and the FDR method was used to adjust for several trials. Finally, we performed a multiple Cox regression analysis to adjust for covariates in order to determine independent survival risk factors.

Results

Autophagy Status and Autophagy-Related DEGs in CRC

The discovery cohort was composed of 574 colorectal cancer patients, and their relapse free survival (RFS) data was retrieved from the GSE39582 data set. Table 1 shows the clinicopathological characteristics of the patients. Tumor staging: 6.4% in stage I, 46.4% in stage II, 35.8% in stage III, and 10.4% in stage IV. We used the expression matrix of 152 hypoxia marker genes in MSigDB version 7.1 to calculate the Euclidean distance between any two patients in the cohort, and used the non-linear dimensionality reduction algorithm t-SNE to compress it into a two-dimensional point. Results identified three patient groups and all patients are each assigned to their closest group (Figure 1A). The patients were classified into three groups: cluster 1(172), cluster 2(189), and cluster(213). Results indicated that there were statistically significant differences in survival rate between these three groups (log rank test, p<0.001). Patients in the first group had the best relapse free survival rate. However, patients in group 2 had the worst prognosis (Figures 1C and D). We further investigated the changes in the expression of genes that regulate autophagy in the two groups of patients. Two gene sets were used, including genes that positively regulate autophagy (120 genes) and genes that negatively regulate autophagy (84 genes). Among the 120 genes associated with “positive regulation of autophagy”, 44 (36.67%) were highly expressed in autophagy^type2 group, and 33 (27.5%) were highly expressed in autophagy^type1 group. Moreover, 30 (35.71%) of the 84 genes associated with “negative regulation of autophagy” were overexpressed in the autophagy^type2 group, and 14 (16.67%) were overexpressed in the autophagy^type1 group. The results above suggest that the identified group has a significant correlation with autophagy status. To obtain autophagy-related DEGs, we compared the expression profile of autophagy^type1 group and autophagy^type2 group. Our results identified a total of 39 DEGs (Table S1) associated with autophagy (Figure 1B). Among them, 26 (67.7%) were overexpressed in the autophagy^type2 group, and the survival rate of patients in the group was comparatively poor. We speculate that these genes are risk factors associated with autophagy. In the autophagy^type1 group, 13 DEGs were highly expressed, and the prognosis of patients was good, indicating that these genes are protective components associated with autophagy. In general, most of the DEGs associated with autophagy are considered to be risk factors.

Table 1 Baseline Characteristics of Patients in the Discovery Cohort

Figure 1 Identification of autophagy and stroma status and autophagy-stroma-related DEGs. (A) Dot plot for three distinct clusters identified by t-SNE algorithm based on 152 autophagy hallmark genes. (B) Heatmap showing expression profiles for autophagy-related DEGs with comparison between type1 and type2 groups. (C) Kaplan–Meier plot of overall survival for patients in three clusters. (D) Kaplan–Meier plot of overall survival for patients in type1 and type2 clusters. (E) Kaplan–Meier plot of overall survival for patients in high-stroma and low-stroma groups. (F) Heatmap showing expression profiles for stroma-related DEGs with comparison between high-stroma and low-stroma groups.

Stroma Status and Stroma-Related DEGs in CRC

We classified the 574 patients according to the condition of stromal cells in tumor tissues and used the “xCell” R package to calculate stromal cell score. We then divided the patients into two groups (the high stromal cell score group and the low stromal cell score group) based on the median of the stromal cell score. Results showed that the prognosis of patients in the stroma^low group was better than that in the stroma^high group (Figure 1E). In addition, we compared the stroma^low group and stroma^high group in order to determine stroma-related DEGs. In total, 19 stroma-related DEGs (Table S2) were screened (Figure 1F), and all the 19 genes were over-expressed in the stroma^high group. Moreover, the patients in the group had relatively poor survival rates. Therefore, we speculate that these DEGs are immune-related dangerous DEGs. In general, most of the DEGs associated with high-stroma cells are considered risk factors.

Autophagy-Stroma-Related Prognostic DEGs in CRC

We further combined the above autophagy and stroma status into a two-dimensional index, and divided the patients into three groups: autophagy^type1/stroma^low group, autophagy^type2/stroma^high group, and the mixed group. Survival analysis showed that there were significant differences among the three groups of patients, where autophagy^type1/stroma^low group had the best prognosis, while the autophagy^type2/stroma^high group had the worst survival rate (Figure 2A). We further compared the expression profiles of autophagy^type1/stroma^low and autophagy^type2/stroma^high groups in order to obtain autophagy-stroma related DEGs. In total, 140 autophagy-stroma related DEGs (Table S3) were identified, of which 25 were overexpressed in the autophagy^type1/stroma^low group. Patients in the group had better survival rates and were defined as autophagy-stroma related protective DEGs. On the other hand, 115 DEGs were overexpressed in autophagy^type2/stroma^high with poor patient prognosis, and thus they were defined as risk DEGs. In addition, we performed Gene Set Enrichment Analysis (GSEA) on the differential genes (Figure 2B), and the results showed that the autophagy^type2/stroma^high group was highly activated in the EMT pathway and inhibited a subgroup of genes regulated by MYC, genes encoding cell cycle related targets of E2F transcription factors, and oxidative phosphorylation pathway. We also overlapped the autophagy (or stroma) related DEGs and the autophagy-stroma related DEGs obtained above to establish a risk gene set (the protective gene set did not have genes that met the requirements). Finally, a total of nine risk genes were identified (Figure 2C). Furthermore, we performed univariate Cox regression analysis on all key risk genes. Among the nine risk factors, six genes had a significant effect on the prognosis of patients.

Figure 2 Identification and biological function of autophagy-stroma-related protective and risk DEGs. (A) Kaplan–Meier plot of overall survival for patients in three groups by combining the autophagy and stroma status. (B) Representative Gene Ontology terms enriched by the autophagy-stroma-related protective and risk DEGs. P-values were adjusted by false discovery rate. (C) Prognostic values of 6 hub genes in Colorectal cancer. Kaplan–Meier analyses were performed based on the median levels of hub genes.

Autophagy-Stroma-Based Gene Classifier and Prognosis Classifier in CRC

We used the LASSO-Cox regression model to determine the most ideal prognostic gene markers from nine autophagy-stroma-related prognostic DEGs that were obtained above (Figures 3A and B). Results indicated that the optimal gene signature was composed of three genes (TAGLN, TNS1, and SFRP4) and the corresponding coefficients were calculated. It is worth noting that the three characteristic genes are all risk genes. The risk score of each CRC patient was calculated as follows: Risk Score=0.01942066×TAGLN expression + 0.007341004×TNS1 expression + 0.03543974×SFRP4 expression. The risk curve and characteristic DEGs expression pattern, based on the risk score of these three genes, are shown in Figure 3D. Comparison of the survival rate showed that the prognosis of patients in the high-risk group was significantly worse than that in the low-risk group (Log rank test P<0.001) (Figure 3C)

Figure 3 Autophagy-stroma-based prognosis classifier. (A) LASSO coefficient profiles of the autophagy-stroma-related prognostic DEGs. (B) Three-fold cross-validation for tuning parameter selection in the LASSO model. The partial likelihood deviance is plotted against log (λ), where λ is the tuning parameter. Partial likelihood deviance values are shown, with error bars representing SE. The dotted vertical lines are drawn at the optimal values by minimum criteria and 1-SE criteria. (C) Survival analysis of high-risk group and low-risk group in the training set. (D) Plot depicting the calibration of the nomogram in terms of the agreement between predicted and observed outcomes. Nomogram performance is shown by the plot relative to the dotted line, which represents perfect prediction. (E) Validation of the model in different data sets.

Autophagy-Stroma-Based Prognosis Signature and Clinicopathologic Characteristics in CRC

Table 1 shows patients’ clinicopathological characteristics in the low-risk group and the group of high risk. Results showed that the gender and age of the two groups were comparable. The group of high risk had more patients with advanced tumors, while the group of low risk had more patients with early tumors (p=0.0119, chi-square test). Moreover, the high-risk group tended to have type 2 autophagy status (p<0.001, χ² test) and patients with high-stroma status (p<0.001, χ² test,), which is consistent with results described above. Univariate analysis in the cohort showed that autophagy and stromal cell status, and prognostic classification were important factors related to survival (p=0.0065, p=0.0081, and p<0.0001, respectively, Log rank test). On the other hand, prognostic classification was identified as an independent prognostic factor by multivariate analysis, independent of tumor stage (Table 2). This suggests the role of stroma cells and autophagy in the prognostic stratification of colorectal cancer.

Table 2 Multivariate Cox Regression Analyses of Risk Factors for Relapse Free Survival

Validation of Autophagy-Stroma-Based Prognosis Classifier

The prognostic classification of autophagy-stroma was further verified using four independent cohorts, including GSE17538, TCGA, GSE38832, and the CMU cohort. The clinicopathological characteristics of the patients are listed in Table 3. The prognostic classifier consistently divided the three validation cohorts into a high-risk group with a poor survival rate and a low-risk group with a good prognosis (GSE17538 group p=0.0052, GSE38832 group p=0.028, TCGA group p=0.046, and CMU group p=0.027, Log rank test) (Figure 3E), which was consistent with results obtained in the discovery cohort. In addition, survival analysis of the four independent cohorts and a verification cohort composed of these four cohorts showed that there was significant difference between the high-risk group and the low-risk group (p<0.0001). Multivariate analysis showed that prognostic classification is an independent prognostic factor in the validation cohort (p=0.00256) and was parallel to the training set (Table 2).

Table 3 TNM Characteristics of Patients in Discovery Cohort and Three Independent Validation Cohorts

The Risk Score Was Significantly Higher in Normal Tissues Than in the Cancer Group, but the Score Gradually Increased with Tumor Development

Furthermore, we analyzed the expression levels of the three key genes of the model constructed using the TCGA pan-cancer data, and found that TNS1 and TAGLN were significantly higher in normal tissues than in tumor tissues (Figures 4A–C). We then compared the scores of 40 normal tissues in the TCGA data set (70 cases of TNM1 stage, 166 cases of TNM2 stage, 107 cases of TNM3 stage, and 65 cases of TNM4 stage, making a total of 448 samples) (Figure 4D). Results showed that the scores in the normal group were significantly higher than those in the cancer group, but the score increases gradually as the cancer develops (increased TNM staging), thereby indicating a worse prognosis.

Figure 4 Differential expression of hub genes in colorectal cancer tissues and para-cancerous tissues and the difference of model scores in different TNM stages of colorectal cancer. (A) Differential expression of TNS1 in colorectal cancer tissues and para-cancerous tissues. (B) Differential expression of TAGLN in colorectal cancer tissues and para-cancerous tissues. (C) Differential expression of SFRP4 in colorectal cancer tissues and para-cancerous tissues. (D) Difference of model scores in different TNM stages of colorectal cancer.

Discussion

In recent years, researchers have focused on the tumor microenvironment. In fact, the tumor microenvironment has been used to construct in vivo and in vitro models for colon cancer research,²⁸ but the clinical application is less. Only a few drugs, such as cetuximab and Avastin, have been used in clinic.⁶ A previous study reported that stromal cells play an important role in the recurrence and metastasis of colorectal cancer.²⁹ Many studies have also confirmed that stromal cells can promote the recurrence and metastasis of colorectal cancer. For example, Hu et al³⁰ reported that stromal cell secreted exosomes promote metastasis and chemotherapy resistance by enhancing cell stemness and epithelial-mesenchymal transition in colorectal cancer. In addition, Ouahoud et al³¹ revealed that bidirectional tumor/stroma crosstalk promotes metastasis in mesenchymal colorectal cancer. Several studies have also shown that stromal cells help to promote cancer cell survival and can protect them against anticancer therapy.^32–35 In this study, results showed that patients with high stromal cell score were more likely to relapse.

During tumor development, the autophagic process appears to function as a tumor suppressor and limits tumorigenesis.¹⁷ The study by Yue et al³⁶ reported that beclin 1 is a critical component of mammalian autophagy and it establishes a role for autophagy in tumor suppression. On the other hand, Takamura et al³⁷ found that autophagy is important for the suppression of spontaneous tumorigenesis through a cell-intrinsic mechanism, particularly in the liver, and that p62 accumulation contributes to tumor progression. Down-regulation of ATG6, a crucial autophagy-related gene, is associated with susceptibility to thyroid cancer³⁸ and colorectal cancer.³⁹ Collectively, these results indicate that autophagy is inhibited during tumorigenesis. There is also a corresponding decrease in model score during the transformation of colorectal cancer from carcinoma in situ to stage 1. However, the enhancement of autophagy may cause adverse effects during the process of cancer progression. For example, Rangwala et al⁴⁰ has shown that inhibiting autophagy can treat advanced solid tumors and melanomas. Vogl et al⁴¹ reported that inhibition of autophagy can treat relapsed/refractory myeloma. Moreover, Cui et al⁴² revealed that radiation induces autophagic cell death in breast cancer cells via the p53/DRAM signaling pathway. These results suggest that autophagy has a controversial effect on cancer. Our results indicated that the model score gradually increases as the cancer progresses, and the patient’s prognosis worsens. Previous studies have reported that cancer cells induce and exploit elevated levels of autophagy in stromal cells for their aberrant growth.^43–45 Therefore, we speculate that the interaction between autophagy and stroma cells can provide prognostic value for patients with colorectal cancer. In this study, we combined autophagy and stromal cell status into the existing clinicopathological characteristics and staging system using a series of systematic analysis, and developed a new genetic marker, with the overarching goal of improving the prognosis of colorectal cancer.

Firstly, we identified genes that are mostly associated with colorectal cancer recurrence under different autophagy states. We classified 574 patients in the GSE39582 data set into three groups using t-SNE, according to the autophagy marker gene set in the MSigDB database. Survival analysis showed that there were significant differences in prognosis among the three groups. Secondly, we chose the two groups with the largest difference, and used them to calculate the correlation between the two groups of patients and the gene set that regulates autophagy. Furthermore, we verified the difference between the two groups. Finally, the genes between the two groups were analyzed to confirm the DEGs associated with autophagy.

Subsequently, we used xCell to score the stroma status of 574 patients, and divided the patients into two groups (high-stroma and low-stroma) according to the score. Survival analysis of the two groups found that the prognosis between the two groups was significantly different. Finally, the genes between the two groups were analyzed to confirm DEGs associated with stroma status.

Next, we divided the patients into autophagy^type1/stroma^low group, autophagy^type2/stroma^high group, and the mixed group. Survival analysis showed that there were significant survival differences among the three groups. We then performed difference analysis between the autophagy^type1/stroma^low group and autophagy^type2/stroma^high group in order to obtain autophagy-stroma-related DEGs.

After obtaining the autophagy/stroma and autophagy-stroma related DEGs, the intersection of the three groups was chosen as the hub genes and used to build the model. In total, we obtained nine genes (SFRP4, ACTG2, SYNPO2, MAB21L2, MGP, SFRP2, CNN1, TNS1, and TAGLN), which are risk genes indicating a poor prognosis when highly expressed. For further survival analysis, six genes with p<0.05 were selected (TAGLN, MGP, TNS1, CNN1, SFRP4, and MAB21L2) and introduced into the Cox model with “glmnet” R package, risk score=0.01942066×TAGLN expression + 0.007341004×TNS1 expression + 0.03543974×SFRP4 expression. Moreover, the 574 patients with colorectal cancer in the GSE39582 data set were scored, and the “maxstat” R package used to determine the optimal cut-point for one or multiple continuous variables at once. Finally, it was verified in 553 patients using GSE17538 (200), GSE38832 (92), TCGA (179), and CMU (82) data sets.

Interestingly, the established model also had some implications for revealing the pathogenesis of colorectal cancer in addition to predicting the recurrence of colorectal cancer. We found that the expression of model key genes, TNS1 and TAGLN, in normal tissues was significantly higher than that in colorectal cancer tissues. However, high expression of TNS1 and TAGLN in colorectal cancer predicted worse prognosis. Moreover, our model had a significantly higher score in normal colorectal tissues than in the cancer group, but the score increased gradually with tumor progression. Therefore, we have reason to believe that normal colorectal tissue experienced a process of decreased stromal cell expression and autophagy changes during carcinogenesis. However, increased expression of stromal cells and changes in autophagy during tumor progression provided favorable conditions for tumor recurrence and metastasis.

Conclusion

In summary, we have established and verified a feasible prognostic model of colorectal cancer based on the autophagy and stromal cell characteristics of patients. The model can be used to evaluate recurrence risk of cancer patients, and the genes in the model provide potential targets for targeted treatment of colorectal cancer. In addition, the mechanism of the occurrence and development of colorectal cancer were reasonably forecasted.