DOI: https://doi.org/10.21203/rs.3.rs-135628/v1
Background: Epithelial-to-mesenchymal transition (EMT) greatly affects the progression and metastasis of colorectal cancer (CRC) which is one of the most malignant tumors. The goal of this study was to evaluate the prognostic value of EMT-related genes and develop a prognostic nomogram for patients’ overall survival (OS) with CRC.
Methods: This retrospective study analyzed gene expression profiles and clinical data in 1099 samples of CRC patients from the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Risk signature based on public databases was created and examined for its predictive value in CRC patients. The independence and correlation of risk signature with clinicopathological parameters were respectively evaluated through univariate/multivariate analysis and subgroup analysis. A nomogram was established according to risk signature and clinical parameters and validated for its predictive ability. The correlation of risk signature and infiltrating immune cells in CRC was assessed.
Results: A four EMT-related genes (ERGs) based signature was constructed and validated for its ability in prognostic prediction in CRC patients. The results showed the ERGs signature as an independent prognostic factor regarding different clinicopathological parameters including age, gender and clinical stage. Besides, the ERGs signature also performed well in subgroup analysis with T3-T4, N-, M0 and early disease stage. A nomogram based on the ERGs signature combining clinical parameters was constructed which showed a promising predictive ability in prognosis of CRC patients’ OS. In addition, we found a positive correlation of the ERGs signature with six different types of immune cells.
Conclusion: A novel ERGs signature was identified in this study and ERGs signature can be useful in predicting the prognosis of CRC patients.
Colorectal cancer (CRC) is the most common cause of gastrointestinal disease related death, ranking third of modality in malignant tumors(1). With a high metastasis-recurrence rate, CRC presents a dismal prognosis for most patients. According to the latest global cancer statistics in 2018, approximately 1.8 million people are diagnosed with CRC and 900,000 patients die from the disease each year(2). Patients with CRC typically undergo surgical treatment, and advancing surgical methods have contributed to ameliorated overall prognosis(3–5). Besides, with the affiliation of adjuvant chemoradiotherapy the overall survival and relapse-free survival for CRC patients have been improved profoundly(3, 6). Additionally, neoadjuvant chemotherapy can increase the percentage of patients who are eligible for resection although the improvement of OS from neoadjuvant chemotherapy remains controversial(7–9). However, uncontrolled metastasis and recurrence still occur in CRC patients despite radical curative surgery and chemoradiotherapy and account for the main cause for the death of CRC patients(10) and the current approaches of assessing CRC patients’ prognosis are not precise enough. Thus, more effective prediction of prognosis for CRC patients are urgently needed for better clinical management.
Epithelial-to-mesenchymal transition (EMT) is regarded as a crucial regulatory process in mediating invasion and metastasis in cancer. EMT enables epithelial cells to acquire a series of mesenchymal characteristics, like motility and invasiveness(11). And recently, studies have shown that epithelial cells are allowed to gain mesenchymal phenotype in which case cells transiently exhibit both epithelial and mesenchymal features without completely transforming into mesenchymal cells as assumed by the classical definition of EMT(12). In contrast with EMT, mesenchymal-epithelial transition (MET) provides cancer cells the loss of mesenchymal phenotypes and the ability of relocating into a distant metastasis. Increasing evidence presents the importance of controlling the EMT process for preventing progression and metastasis in multiple cancers, including CRC, and methods on treatment of targeting EMT process are also in researching(13). However, the clinicopathological character of EMT in CRC has not been well understood. Therefore, we performed systemic bioinformatic analysis in this study to explore the potential role of EMT in prediction of patient prognosis and discover the possible novel modality targeting EMT process to improve existing treatment.
In this study, we filtrated four EMT-related genes (ERGs) from differentially expressed genes (DEGs) by analyzing gene datasets in Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) data sets. We evaluated the correlation of ERGs signature with clinical characteristics and tumor associated immune cells through which the risk signature may reflect CRC microenvironment and immune response. Moreover, we integrated the ERGs signature with clinical factors to build a prognostic nomogram, which allowed ameliorated prognosis assessment of CRC patients.
Study design of this research is presented in the flow chart (Fig. 1A). Level 3 RNA sequencing data (RNA-seq) and corresponding clinical information of 514 CRC patients was obtained from the TCGA database (https://portal.gdc.cancer.gov/repository)(14). The gene expression profiles were normalized by scale method in “limma” R package. RNA-seq data and paired clinical information of another 585 CRC patients was downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/)(15) and corresponding dataset was GSE39582. Normalized count values were used. Both TCGA and GEO database are publicly available, thus the present study was exempted from the approval of local ethics committees. Criteria for study inclusion were: (1) Repeated tumor samples in the same patient were removed. (2) Patients with unknown survival status and follow-up information, and those who died within a follow-up period of 30 days were excluded. (3) Patients with unknown disease stage or grade were excluded. Consequently, 416 patients meeting the criteria were included in the discovery set and 532 patients in validation set. The TCGA discovery set was used to establish an ERG signature to predict patients’ prognosis, while the GEO set were used for validation.
200 ERGs were downloaded from the Gene Set Enrichment Analysis (GSEA) database (http://software.broadinstitute.org/gsea/index.jsp)(16).
Differentially expressed mRNA between tumor and adjacent normal tissue samples were screened in TCGA cohort by the “limma” package of R software 4.0.2. |log2 fold change|>1 and false discovery rate (FDR) < 0.05 were set as the cut-off criteria. DEGs in paired samples from TCGA datasets, characterized by upregulation and downregulation, were intersected with GSEA ERGs set to obtain differentially expressed ERGs.
After identifying hub ERGs, we performed univariate Cox analysis of overall survival (OS) to screen ERGs with prognostic value. Subsequently, we performed least absolute shrinkage and selection operator (LASSO) Cox regression analysis in TCGA discovery set to select the best gene model for predicting prognosis in CRC patients. LASSO Cox regression analysis was performed using R package “glmnet” and the optimal values of penalty parameters were determined by 10-fold cross-validation. Furthermore, Kaplan-Meier (K-M) survival curves were used to assess survival difference between high- and low-risk groups using the “survminer” package in R software. The risk score model for each patient was determined by a linear combination of gene expression weighted by the regression coefficient from LASSO Cox regression analysis. Then, we used the median value of the risk score as cut-off value which was next used to divide patients into high- and low-risk group. To assess the predictive accuracy of the ERG signature time-dependent Receiver Operator characteristic (ROC) analysis was performed with area under the curve (AUC) calculated at different cut-off time using the “survivalROC” package in R software.
The same ERGs signature model was used in GEO set for external validation. The risk score was calculated for each patient in validation set and X-tile software(17) was utilized for determining the optimum cut-off point. Survival analysis was performed in GEO cohort to validate the prognostic ability of ERG-signature. Additionally, time-dependent ROC curves were employed to evaluate the prediction accuracy of ERG-signature model in validation set.
To determine the independence of the prognostic model from other clinical parameters, univariate and multivariate analysis of risk-score with age, gender, disease stage was performed. The relationship between the risk score and clinicopathological parameters was showed in forest plots.
To evaluate whether the prognostic model was associated with other clinicopathological parameters including age, gender, disease stage, pathological T stage, pathological N stage and pathological M stage, patients were divided into subgroups of age ≤ 65 years and age ≥ 65 years, female and male, pathological tumor stage I-II and stage III-IV, pathological T1-T2 and T3‐4, pathological N− and N+ (N1‐3), without metastasis M0 and with metastasis M1. K-M survival analysis was performed for the subgroups of different clinicopathological parameters.
A nomogram integrating the ERG signature and various clinicopathological factors was established for predicting patients’ prognosis using the “rms” package in R. Calibration curves were applied to assess the accuracy of the nomogram.
The fraction data of six immune cells, including B cell, CD4+ T cell, CD8+ T cell, Neutrophil, Macrophage and Dendritic, was downloaded from the Tumor Immune Estimation Resource (TIMER) website (https://cistrome.shinyapps.io/timer/)(18) and their correlation with risk score was calculated.
The heat map was generated using the “pheatmap” package in R. Student’s t-test was used to compare subgroups and paired data. Pearson’s chi-square test was applied to analyze differences between the discovery and validation sets and the association between the risk score and clinical parameters. K-M survival curves were compared using log-rank test. Univariate Cox regression was conducted to estimate the hazard ratio (HR) for different factors. Multivariate Cox regression was performed to determine independent factors. All statistical analyses were conducted by R software. P < 0.05 was considered significant.
Gene expression profiles of CRC patients from TCGA database were extracted and intersected with 200 EMT-related genes downloaded from GSEA dataset. Then, we compared gene expression levels of the filtrated genes between tumor tissues and corresponding normal tissues and identified 62 DEGs of which 23 were down-regulated and 39 were up-regulated. According to the DEGs a heatmap was drawn and the distribution of DEGs expression was displayed in the Volcano plot (Fig. 1B&C).
We used Univariate Cox analysis to identify ERGs with prognostic value and seven ERGs with highest Hazard ratios were screened out and a forest plot was performed (Fig. 1D). To identify ERGs with most valuable prognostic efficacy, we performed the least shrinkage and selection operator (LASSO) Cox regression model in the TCGA discovery set. Consequently, we obtained four ERGs, namely PCOLCE2, OXTR, BDNF, SERPINE1 to construct our gene signature. Besides, all the four ERGs were negatively correlated with patient survival.
Based on the respective expression levels, we constructed a risk score formula: Risk score = 0.246592982*Exp(BDNF) + 0.070329939*Exp(OXTR) + 0.052615049*Exp(SERPINE1) + 0.276502574*Exp(PCOLCE2). (Exp: expression level of the gene)
Expression levels of the four ERGs in TCGA cohorts were exhibited in a heatmap (Fig. 2A). Patients were separated into high- and low-risk score groups according to a median cut-off value of 0.546. The distribution of risk scores and survival time were shown in Fig. 2B&2C. Poorer prognosis was presented in the high-risk group than the low-risk one. A time-dependent ROC analysis was performed to show the prognostic accuracy of the four-gene signature. The AUC was 0.639 at 1 year, 0.664 at 3 year and 0.603 at 5 year (Fig. 2D). K-M survival analysis corroborated that patients in the high-risk group showed a poorer overall survival rate than the low-risk(p < 0.01) (Fig. 2E).
To further verify the prognostic efficacy of four-ERG signature, we examined the prognostic potential in a GEO validation set (Fig. 3A). Risk-score of each patient was calculated by the same formula but the patients in GEO cohort was separated into high- and low-risk groups using an optimum cut-off point (2.6) determined by X-tile software (Fig. 3B). Consistent with the result in discovery set, high-risk patients in validation set had a poorer prognosis than low-risk group (Fig. 3C). Besides, ROC curve performed well in the validation set with AUC of 0.615 at 1 year, 0.583 at 3 year and 0.574 at 5 year (Fig. 3D). Similarly, a lower survival rate was exhibited in high-risk patients from validation set by survival analysis (Fig. 3E).
To verify four-ERG signature was an independent prognostic factor regarding clinical characteristics, we performed univariate and multivariate cox analysis (Fig. 4). The result showed that the risk score could predict the prognosis independently regarding to age, gender and clinical stage.
To further evaluate the prognostic value of four-ERG signature in specific clinical conditions, subgroup analysis between risk score and different clinical parameters was performed. The result suggested that risk signature stayed powerful in either younger or older, male or female. However, the result was quite intriguing in subgroups of tumor node metastasis (TNM) stage and disease stage. Our risk signature showed a convincing predictive power in subgroup of T3-T4 (p = 0.001), N− (p < 0.001), M0 (p < 0.001) and early disease stage (stage I-II, p < 0.001) while high-risk group also tended to have a poorer prognosis but the difference didn’t reach significance in T1-2 (p = 0.282, N+ (p = 0.213) and advanced disease stage (stage III-IV, p = 0.290). Besides, the risk score showed no efficient predictive ability in the metastasis (M1) subgroup (p = 0.714) (Fig. 5).
With the aim to develop a feasible tool for clinical predicting OS of CRC patients, we constructed a nomogram with the risk score (based on the ERGs signature), age, gender and disease stage in the TCGA dataset (Fig. 6A). Calibration plots showed agreeable concordance between observed and predicted rates of 1-, 3- and 5-year OS. Besides, time-dependent ROC curves were used for assessing the prognostic predictive ability of the nomogram and other predictors. The result indicated that the nomogram owned satisfying accuracy of 1-, 3- and 5-year OS in combination with age, gender and disease stage (Fig. 6B).
To assess the efficacy of ERGs signature reflecting the tumor immune microenvironment status, we evaluated the risk score and the representation of different immune cells. The result showed that ERGs signature was positively associated with CD4+ T cells (Cor = 0.459, p < 0.01), CD8+ T cells (Cor = 0.247, p < 0.01), B cells (Cor = 0.223, p < 0.01), dendritic cells (Cor = 0.453, p < 0.01), macrophages (Cor = 0.513, p < 0.01) and neutrophils (Cor = 0.379, p < 0.01) (Fig. 7). These data suggested that the four-ERG signature could serve as an indicator of immune status of CRC.
In consideration of the importance of epithelial-mesenchymal transition in malignancy metastasis(19), developing meaningful gene signatures to monitor the EMT status and motility of tumor cells in cancer patients is significant not only to identify reliable prognostic biomarkers but, if applied properly, to stratify patients at high-risk of recurrence, metastasis or chemoresistance who might benefit from additional therapy. In the present study, we validated a prognostic signature, based on four ERGs, which proved to be a reliable indicator of EMT status and could identify CRC patients with unfavorable prognosis. Moreover, our prognostic signature can further stratify patients sharing specific clinicopathological factors (e.g., age, gender and disease stage) into subgroups with different survival outcomes. Integrating these findings, we constructed a nomogram, incorporating the ERGs signature’s risk score and clinical characteristics, which showed good performance for predicting survival in patients with CRC.
Although several gene signatures representative of tumor EMT status have been proposed with potential clinical applicability in several cancers(16, 20), the impact of EMT status in CRC progression and prognosis prediction remains to be fully explored. Subgroup analysis combining our ERGs signature risk score and clinicopathological parameters revealed that a higher risk score was more strongly correlated with poorer prognosis. This is owing to high risk score reflecting the possibilities of enhancing the motility and invasiveness of tumor cells contributing to tumor progression and metastasis. Noteworthy, we also confirmed that the ERG signature remained an independent prognostic factor with multivariate analysis after adjusting for clinicopathological variables.
Cancer progression and metastasis is often linked to altered plasticity of tumor cells(21, 22), while EMT plays an unignorable role in this process(23). Accordingly, aberrant expression of transcription factors, cytokines, epithelial surface markers and mesenchymal features has been shown to contribute to the migration of CRC cells(24). According to previous studies, the four ERGs in our novel CRC signature encode proteins with active participation in angiogenesis, developing resistance and chemotaxis, increasing cancer cell growth, proliferation, survival, migration and epithelial to mesenchymal transition, and decreasing anoikis, relapse and chemotherapeutic sensitivity.
Brain-derived neurotrophic factor (BDNF) is highly expressed in brain as member of neurotrophins which were initially known to function as growth factors in nervous system(25). However, recent studies have discovered that BDNF and its receptor tropomyosin receptor kinase B (TrkB) were upregulated in a wide range of cancers including CRC(26). Growing evidence showed that BDNF/TrkB pathway was closely involved in CRC proliferation, invasion, migration and poor prognosis(27–29). Besides, positive results render BDNF/TrkB pathway to be a promising therapeutic target for CRC treatment(27, 30).
Oxytocin (OT) is viewed as a hypothalamic neuropeptide and classically functions in uterine contraction and milk ejection(31). Interestingly, some findings suggested that the gene encoding OT receptor (OXTR) also expressed in several types of cancer cells including breast, ovarian and prostate cancer regulating proliferation(32, 33). The relation between OXTR and CRC hasn’t been validated though it was reported that oxytocin receptors were expressed throughout gastrointestinal tract(34).
Plasminogen activator inhibitor 1 & 2 (PAI-1 & PAI-2), also known as serpin family E member 1 & 2 (SERPINE1 & SERPINB2), are the main inhibitor of urokinase plasminogen activator (uPA) and tissue plasminogen activator (tPA) which constitute plasminogen activator (PA) system. PA system in turn acts crucially in extracellular matrix (ECM) remodeling and further plays an important role in tumor progression(35). However, SERPINE1 was found associating with tumor progression and metastasis and identified as an indicator of poor prognosis in a series of cancers(36–38). Similarly, the character of SERPINE1 (or PAI-1) in colorectal cancer was also identified, as a prognostic biomarker(39, 40).
PCOLCE2 encodes procollagen C-proteinase enhancer 2 which is collagen-binding protein capable of binding at multiple sites on the triple helical portions of fibrillar collagens and competing for such binding with procollagen C-proteinases(41). In consistent with our research, PCOLCE2 has already been identified as a diagnostic and prognostic biomarker in CRC patients(42). However, the underlying mechanism of PCOLCE2 engaging in the progression and metastasis of CRC remains to be elucidated and we for the first time propose the idea that PCOLCE2 is associated with EMT process of CRC cancer cells.
Accumulating evidence shows that tumor cells undergoing EMT or acquiring mesenchymal phenotype might have contribution to the component of tumor infiltrating immune cells (TIICs) shifting to an immunosuppressive pattern resulting in immune escape and immune resistance(43, 44). Studies have also showed cross-talk between EMT and different immune cells render the availability to predict the presentation of TIICs by assessing EMT status(20, 45–47). In this study, ERG signature showed a positive correlation with immune cells including CD4+ T cells, dendritic cells and macrophages. In the light of our knowledge, CD4+ T cells and macrophages contribute greatly to the formation of immunosuppressive microenvironment, however, different subpopulations of these immune cells exhibit diverse functions and prognostic roles in CRC, thus immune cells can’t be simply regarded as risk or protective factors to CRC(46). Nevertheless, our study showed a certain change in infiltrating immune cells regarding to the ERGs, indicating the potential value of further investigating the relationship between the ERGs and subpopulations of different tumor-associated immune cells and their prognostic and therapeutic role in CRC.
In summary, we identified a gene signature for CRC based on four differentially expressed ERGs, which could act as an indicator for tumor EMT status and to classify patients into groups with distinctly different prognosis. Additionally, we constructed a composite prognostic nomogram for prediction of CRC patients’ prognosis integrating ERGs signature and clinical factors, which may serve as a promising clinical tool. Limitation of our study lay in its retrospective nature and the need to define the molecular mechanisms by which the ERGs in our signature impact CRC progression and metastasis. Nevertheless, besides aiding current prognostic efforts, we hope the data provided will lead to further studies aimed at exploring novel therapeutic strategies for CRC.
CRC: colorectal cancer
DEGs: differentially expressed genes
EMT: epithelial-to-mesenchymal transition
ERGs: EMT-related genes
OS: overall survival
TCGA: The Cancer Genome Atlas
GEO: Gene Expression Omnibus
MET: mesenchymal-epithelial transition
GSEA: Gene Set Enrichment Analysis
FDR: false discovery rate
LASSO: least absolute shrinkage and selection operator
K-M: Kaplan-Meier
ROC: time-dependent Receiver Operator characteristic
AUC: area under the curve
HR: hazard ratio
Cor: correlation coefficient
BDNF: brain-derived neurotrophic factor
TrkB: tropomyosin receptor kinase B
OT: oxytocin
OXTR: OT receptor
PAI-1 & PAI-2: plasminogen activator inhibitor 1 & plasminogen activator inhibitor 2
SERPINE1 & SERPINB2: serpin family E member 1 & serpin family E member 2
uPA: urokinase plasminogen activator
ECM: extracellular matrix
PCOLCE2: procollagen C-proteinase enhancer 2
TIICs: tumor infiltrating immune cells
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Competing interests
The authors declare that they have no competing interests
Funding
This work was supported by grants from the National Natural Science Foundation of China (Grant No. 81372553)
Authors' contributions
SG and CY made substantial contributions to conception and design of the research. SH, CY and FC carried out data collection, analysis. SH, CY and FC provided interpretation of the results. SH and CY wrote the paper. SH, CY, FC, XH, and YZ edited the manuscript and provided critical comments. All authors read and approved the final manuscript.
Acknowledgements
Not applicable