1. Cancer Genome Atlas (TCGA) analysis revealed metastasis-related differentially expressed genes
We downloaded raw mRNA expression data (FPKM) of processed COAD from the TCGA database and obtained the metastasis through clinical indicator ( M = 1: metastasis; M = 0: no metastasis ). Differential expression analysis was performed between colon cancer metastasis patients and colon cancer non-metastasis patients through the 'Limma' package. From the research results, we can see that a total of 553 genes that might be related to metastasis were differentially expressed ( | log2FC | ≥ 1 and p < 0.05 ). Among them, 331 genes were up-regulated and 222 genes were down-regulated( Fig. 2.A ). Then 552 differentially expressed genes were screened by gene expression value ( genes whose expression value was not 0 or the average expression value was more significant than 0.3 in 50% samples ). GO and KEGG enrichment analysis showed that These differential genes are enriched in GO pathways such as antigen processing presentation (antigen processing presentation), tumor necrosis factor response (response to tumor necrosis factor), and MHC protein complex binding (MHC protein complex binding)( Fig. 2.B ). In the process of KEGG enrichment, there are Antigen processing and presentation ( antigen processing and presentation ), Primary immunodeficiency ( primary immune deficiency ), TNF signaling pathways, and a large number of genes are enriched in metabolic-related pathways ( Fig. 2.C ).
Moreover, the basic characteristics of the Training and Test Sets are listed in Table 1.
Table 1
The Characteristics of Patients in the Training and Test Sets
|
TCGA testing
|
TCGA training
|
p
|
n
|
80
|
318
|
|
futime (mean (SD))
|
780.88 (890.88)
|
710.01 (685.21)
|
0.439
|
fustat = 1 (%)
|
12 (15.0)
|
58 (18.2)
|
0.606
|
age (mean (SD))
|
68.90 (11.89)
|
66.74 (12.58)
|
0.166
|
gender = MALE (%)
|
35 (43.8)
|
176 (55.3)
|
0.083
|
stage (%)
|
|
|
0.433
|
Stage I
|
12 (15.0)
|
61 (19.2)
|
Stage II
|
3 ( 3.8)
|
23 ( 7.2)
|
Stage IIA
|
29 (36.2)
|
95 (29.9)
|
Stage IIB
|
1 ( 1.2)
|
8 ( 2.5)
|
Stage IIC
|
0 ( 0.0)
|
1 ( 0.3)
|
Stage III
|
1 ( 1.2)
|
13 ( 4.1)
|
Stage IIIA
|
1 ( 1.2)
|
5 ( 1.6)
|
Stage IIIB
|
9 (11.2)
|
42 (13.2)
|
Stage IIIC
|
11 (13.8)
|
18 ( 5.7)
|
Stage IV
|
10 (12.5)
|
35 (11.0)
|
Stage IVA
|
2 ( 2.5)
|
13 ( 4.1)
|
Stage IVB
|
0 ( 0.0)
|
2 ( 0.6)
|
unknow
|
1 ( 1.2)
|
2 ( 0.6)
|
T (%)
|
|
|
0.367
|
T1
|
2 ( 2.5)
|
7 ( 2.2)
|
T2
|
11 (13.8)
|
61 (19.2)
|
T3
|
62 (77.5)
|
209 (65.7)
|
T4
|
3 ( 3.8)
|
24 ( 7.5)
|
T4a
|
2 ( 2.5)
|
11 ( 3.5)
|
T4b
|
0 ( 0.0)
|
6 ( 1.9)
|
M (%)
|
|
|
0.723
|
M0
|
68 (85.0)
|
268 (84.3)
|
M1
|
11 (13.8)
|
39 (12.3)
|
M1a
|
1 ( 1.2)
|
8 ( 2.5)
|
M1b
|
0 ( 0.0)
|
3 ( 0.9)
|
N (%)
|
|
|
0.143
|
N0
|
46 (57.5)
|
198 (62.3)
|
N1
|
13 (16.2)
|
50 (15.7)
|
N1a
|
2 ( 2.5)
|
10 ( 3.1)
|
N1b
|
0 ( 0.0)
|
11 ( 3.5)
|
N1c
|
0 ( 0.0)
|
2 ( 0.6)
|
N2
|
18 (22.5)
|
36 (11.3)
|
N2a
|
0 ( 0.0)
|
5 ( 1.6)
|
N2b
|
1 ( 1.2)
|
6 ( 1.9)
|
2. Construction of the prognosis model for metastatic colorectal cancer patients
In order to make the study more valuable and to explore the key genes of differential gene concentration related to metastasis of colorectal cancer, the relevant data of TCGA-COAD patients were integrated in this study and Cox univariate regression and Lasso regression feature selection algorithm were used to screen feature genes related to colon cancer metastasis ( Fig. 3.A-C ). The results showed that 26 prognostic genes were screened by Cox univariate regression, which was GAL, UCHL1, TRIP10, SERPINE1, SNAI1, BCL10, GSR, PHF2, DNAJB2, LRRC8A, CST6, JAG2, ASAH1, C4orf19, MOGS, GDI1, SNCG, ASRGL1, LEPROTL1, FDFT1, CNOT7, TSC22D3, TNK2, RNASET2, CPT2, PGM2 ( Fig. 3.A ). After the screening of characteristic variables by Lasso regression, 14 selected metastasis related genes were finally determined ( Fig. 3.B-C ). TCGA patients were randomly assigned, and the ratio of training set to validation set was 4 to 1 to construct and verify the Cox proportional hazards regression model. The optimal MR-risk score for each sample was obtained by regression model for subsequent analysis. (MR-Risk Score = ASRGL1 × (-0.200486209) + GSR × (-0.107200303) + ASAH1 × (-0.086956773) + BCL10 × (-0.065018786) + SNAI1 × 0.001969402 + TRIP10 × 0.02686487 + TSC22D3 × 0.063610852 + LRRC8A × 0.072751154 + PHF2 × 0.075148411 + SERPINE1 × 0.077260841 + RNASET2 × 0.110933675 + DNAJB2 × 0.13931055 + UCHL1 × 0.205894794 + GAL × 0.209549353 )( Fig. 3.D ). COX coefficient is shown in Table 2.
Table 2
gene
|
coef
|
hr
|
low.ci
|
upp.ci
|
ASRGL1
|
-0.20049
|
0.559465
|
0.370334
|
0.845184
|
GSR
|
-0.1072
|
0.602317
|
0.445522
|
0.814294
|
ASAH1
|
-0.08696
|
0.650391
|
0.487754
|
0.867258
|
BCL10
|
-0.06502
|
0.586685
|
0.433268
|
0.794424
|
SNAI1
|
0.001969
|
1.424898
|
1.179443
|
1.721434
|
TRIP10
|
0.026865
|
1.548546
|
1.243982
|
1.927677
|
TSC22D3
|
0.063611
|
1.253047
|
1.06009
|
1.481127
|
LRRC8A
|
0.072751
|
1.298438
|
1.099505
|
1.533362
|
PHF2
|
0.075148
|
1.418221
|
1.140282
|
1.763906
|
SERPINE1
|
0.077261
|
1.222669
|
1.100739
|
1.358106
|
RNASET2
|
0.110934
|
1.268938
|
1.062059
|
1.516116
|
DNAJB2
|
0.139311
|
1.42694
|
1.141169
|
1.784273
|
UCHL1
|
0.205895
|
1.380657
|
1.195374
|
1.594658
|
GAL
|
0.209549
|
1.525341
|
1.270088
|
1.831894
|
3. Kaplan-Meier curve and time-varying ROC analysis were performed based on the risk values of training and test data sets
According to the median MR-Risk score (training set risk score median: -0.0374618011019423. The median MR-risk score in the test set was 0.0223075396080433 ). In the TCGA-COAD dataset, patients were divided into two groups, namely high-risk and low-risk groups based on median MR-risk score in the training (Fig. 4.A) and test set (Fig. 4.D ), and studied by Kaplan-Meier curve. We can see that the low risk group has a significantly higher OS than the high risk group ( Fig. 4.B & E ). The results of ROC curve show that the C-index index of training set and test set is 0.72 and 0.76 ( Fig. 4.C & F ), respectively, suggesting that the model has good verification efficiency.
4. External validation of the performance of the MR-risk score model in GEO dataset
In order to ensure the stability of the prediction model, the RNA-Seq data (GSE38832) of COAD patients with processed survival data was retrieved from the GEO database, an external data set, and based on this model, the clinical classification of COAD patients was predicted. In the GSE38832 dataset, patients were divided into two groups, namely high-risk and low-risk groups based on median MR-risk score (Fig. 5.A & B). Kaplan-Meier method was used to evaluate survival differences between the two groups. The results showed that the high risk group in the GEO external validation set had a significantly lower OS than the low risk group( Fig. 5.C ).To verify the accuracy of the model, we used the external data set to analyze the ROC curve of the model, and the results showed that the model had a strong predictive effect on the prognosis of patients( GSE38832 C-index = 0.68 ) ( Fig. 5.D ). At the same time, the expression levels of 14 genes in the prediction model were detected in the GEO external validation set, and it was found that the expression levels of 14 genes were highly consistent in the GSE38832 data set ( Fig. 5. E ) and the TCGA-COAD dataset ( Fig. 5. F ).
5. 14 metastasis-related prognostic features associated with overall survival in patients with colorectal cancer
In order to prove the application value of this prediction model in clinical practice, We first used univariate and multivariate Cox risk regression models to evaluate the impact of MR-risk score on patients' survival( Fig. 6.A-B ). The results showed that MR-risk score was an independent prognostic factor in patients with COAD. Clinically, we often classify the progression and malignant degree of tumors by staging classification so as to achieve the purpose of different surgical and drug treatments for tumors with further advancement and malignant degrees. According to the above, we divide each patient's MR-risk score value evaluated by the prediction model into different clinical commonly used staging classification groups and display the results of each clinical index group in the form of a box diagram ( Fig. 6.C ) and find that the distribution of MR-risk score value in stage, T, M, N and other clinical ( p < 0.05 ) by rank-sum test ( Kruskal. test ). It shows that the MR-risk score obtained by our modeling analysis has good applicability for the grouping of samples, and with the increase of MR-risk score value, it indicates that the clinical tumor staging classification of patients shows a trend of deterioration and the possibility of tumor metastasis increases.In order to play the practical role of the prognostic model, we classified the samples into high risk group and low risk group based on the median of MR-risk score during the study.The results of regression analysis are shown in the form of a line chart( Fig. 6.D ). In order to observe whether the prediction probability of the model is close to the empirical probability, we build the calibration plot by using the predicted five-year survival and seven-year survival ( OS ) of colon cancer( Fig. 6.E ). We also calculated the proportion of cancer metastasis in high versus low-risk score subgroups, results shows that in high-risk subgroup patients have a higher probability of cancer metastasis( Fig. 6.F ).
6. The immune infiltration performance of the model in the high-risk and low-risk Groups
The microenvironment of tumor formation is generally composed of a variety of complex components such as fibrocells and immune cells related to tumor, the extracellular matrix, a variety of growth factors that promote cell growth, cancer cells themselves and related inflammatory factors, and unique physical and chemical components that promote the formation of cancer cells. Tumor microenvironment has a certain or significant influence on the judgment of tumor and patient survival as well as clinical treatment. Studying the association between MR-Risk score and tumor immunoinvasion can further optimize the mechanism effect of MR-Risk score in predicting potential colon cancer molecules.The results showed that the MR-risk score was positively correlated with Macrophages M0, Treg, and B cell naive, and negatively correlated with T cell CD4 memory resting, Eosinophils, Neutrophils and T cell follicular helper ( Fig. 7.A-C). The correlation coefficient scatter plot of T cells gamma delta were not shown because the cell infiltration rate is too low to show the linear correlation between high and low MR-risk score. In order to evaluate the possible sensitivity of MR-risk score to immunotherapy, we examined the correlation between MR-risk score and the expression level of immune checkpoint in tumor tissues( Fig. 7.D). Interestingly, we found a negative correlation between MR-risk score and gene expression at most immune checkpoints in cancer tissues.