Identification and Verification of a 17 Immune-Related Gene Pair Prognostic Signature for Colon Cancer

Abstract

Background. Colon cancer (CC) is a malignant tumor with a high incidence and poor prognosis. Accumulating evidence shows that the immune signature plays an important role in the tumorigenesis, progression, and prognosis of CC. Our study is aimed at establishing a novel robust immune-related gene pair signature for predicting the prognosis of CC. Methods. Gene expression profiles and corresponding clinical information are obtained from two public data sets: The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO, GSE39582). We screened out immune-related gene pairs (IRGPs) associated with prognosis in the discovery cohort. Lasso-Cox proportional hazard regression was used to develop the best prognostic signature model. According to this, the patients in the validation cohort were divided into high immune-risk group and low immune-risk group, and the prediction ability of the signature model was verified by survival analysis and independent prognostic analysis. Results. A total of 17 IRGPs composed of 26 IRGs were used to construct a prognostic-related risk scoring model. This model accurately predicted the prognosis of CC patients, and the patients in the high immune-risk group indicated poor prognosis in the discovery cohort and validation cohort. Besides, whether in univariate or multivariate analysis, the IRGP signature was an independent prognostic factor. T cell CD4 memory resting in the low-risk group was significantly higher than that in the high-risk group. Functional analysis showed that the biological processes of the low-risk group included “TCA cycle” and “RNA degradation,” while the high-risk group was enriched in the “CAMs” and “focal adhesion” pathways. Conclusion. We have successfully established a signature model composed of 17 IRGPs, which provides a novel idea to predict the prognosis of CC patients.

1. Introduction

The incidence of colorectal cancer ranks third but second in terms of mortality around the world, and colon cancer (CC) accounts for seventy percent of colorectal cancer [1, 2]. At present, the treatment strategy is still based on surgery and supplemented by radiotherapy and chemotherapy. However, the overall five-year survival rate of CC is still only sixty-three percent [1], which is mainly due to the lack of early screening [3, 4] and the formation of resistance to postoperative standard chemotherapy drugs [5–7].

In recent years, targeted therapies such as the application of drugs targeting VEGF or EGFR [8], and immunotherapy [9–12] can be used to treat some patients with advanced CC. It should be noted that patients with the same clinical characteristics, pathological types, and treatment get different prognostic results, indicating that the innate genetic heterogeneity of patients has a great impact on clinical and molecular diversity [13].

Besides, the treatment of colorectal cancer is becoming more and more individualized and differentiated, so CC patients need to establish a clinical prognostic model based on the survival rate of patients.

Immunotherapy is the use of drugs to help a person’s own immune system better recognize and destroy cancer cells, which has played an important role in the treatment of CC in the past few years [14]. CC cells prevent the immune system from attacking them by acting on “checkpoint” on immune cells if patients whose CC cells tested positive for specific gene changes, drugs called “immune checkpoint inhibitors” can be used [15, 16]. For example, the application of PD-1 [14, 15] inhibitors and CTLA-4 [14, 17] inhibitors in the treatment of CC has achieved certain results, so immunotherapy has received more and more attention.

The interaction between CC immune genes and their value in the prognosis of CC remains to be further studied. In the present study, through retrospective analysis of the public dataset, the CC immune-related genes were screened, and on the basis of what the immune-related gene pair model was built, the validity and accuracy of its prognostic features were verified.

2. Material and Methods

2.1. Data Collection and Processing

The transcriptome profiling data was downloaded from The Cancer Genome Atlas (TCGA-COAD, https://portal.gdc.cancer.gov/) treated as the discovery cohort (), and the validation cohort was obtained in the form of the microarray from Gene Expression Omnibus (GSE39582, , https://www.ncbi.nlm.nih.gov/geo/).

The gene expression profiles (GEPs) and the corresponding clinical data of the two datasets were gained and processed, respectively. In addition, the normal samples were removed, and only the tumor samples with complete survival information were retained for further analysis.

For the discovery cohort, the collection of transcriptome data, the conversion of Ensembl IDs, and the extraction of relevant clinical data were all done by performing Strawberry Perl (5.30.11). For the validation cohort, converting the probe matrix into the gene matrix was also completed by Perl.

2.2. Immune-Related Gene Expression Data

A list of 2498 immune genes and their action categories was retrieved from ImmPort (https://www.immport.org/shared/home) visited on March 20, 2020. Referring to the list, we gained the expression data of immune-related genes (IRGs) from the transcription matrix of discovery cohort and validation, respectively.

2.3. Prognosis-Related Immune Pairs

The IRGPs and clinical data were analyzed jointly. In our study, the IRGPs were measured on the platform with high variability, which was determined by the median absolute deviation [7] >0.5.

If the expression level of the first IRG is lower than the second one, the score is 0; otherwise, the score is 1. If more than 80% of the samples in the discovery cohort were score 1 or 0, the samples were filtered out, and the remaining IRGPs were left as the initial candidate prognosis-related immune gene pairs ().

2.4. Construction of IRGP Signature Model

After using Lasso-Cox proportional risk regression (1000 simulation iterations), we got the stable IRGPs to construct the final prognosis model (“glmnet” package, version: 3.0-2).

The model signature was presented as the risk

We found the optimal cut-off value through the ROC curve (“survival ROC” package), that is, the point with the greatest sum of sensitivity and specificity. Samples with a higher risk score than the cut-off value were classified as a high-risk group, and vice versa.

2.5. Validation of IRGP Signature

To verify the accuracy of patients’ prognoses stratified by the IRGP signature, we used survival analysis and independent prognostic analysis to see whether there was statistical significance between the high-risk and low-risk groups.

2.6. Immune Cell Infiltration

CIBERSOFT [18] was adopted to compare the expression of immune cells between the two risk groups, which is a versatile computational method for quantifying cell fractions from bulk tissue GEPs. To quantitatively capture deconvolution confidence, CIBERSOFT calculates several quality control metrics, including a deconvolution value.

TCGA-GEPs were uploaded to the CIBERSOFT portal (http://CIBERSORT.stanford.edu/), and the abundance of 22 kinds of infiltrating immune cells in each sample was calculated, including monocytes, macrophages, B cells, and T cells. The correlation of immune cell infiltration between the two risk groups was visualized by boxplot and radar chart.

2.7. Gene Set Enrichment Analysis (GSEA)

In order to observe the functional pathways related to the differentially expressed genes between the two risk groups, we conducted GSEA on the discovery cohort, and the KEGG (Kyoto Encyclopedia of Genes and Genomes) dataset was retrieved from GSEA datasets (“c2.cp.kegg.v7.1.symbols.gmt”, https://www.gsea-msigdb.org/). Finally, the visualization of GSEA results was shown by performing the “fgsea” package and “ggplot2” package in R software ().

2.8. Statistical Analysis Method

All statistical analysis was completed by R software (version 3.6.3), and data extraction and processing were realized by performing Strawberry Perl (version 5.30.2.1).

The differences among groups were compared by Student’s -tests or Wilcoxon rank-sum tests. The Kaplan-Meier method was used for survival analysis, and the “survival” R-package was used for the log-rank test. Cox proportional hazards regression model was used for univariable and multivariable analyses. For all the analyses, a value of less than 0.05 was considered statistically significant.

3. Results

3.1. Construction and Evaluation of Prognostic IRGP Signature

This is a retrospective study of 1037 patients. The TCGA dataset () and GEO dataset () were used as the discovery cohort and the validation cohort, respectively (Table 1).

In the beginning, we screened 326 immune-related genes and gained 12276 immune-related gene pairs (IRGPs). After removing the less variable IRGPs and analyzing together with clinical data, Lasso-Cox proportional hazard regression was used to define the IRGP signature. Then, a total of 17 IRGPs composed of 26 IRGs were applied to construct the risk scoring model. Finally, we calculated the risk score for each patient in the discovery cohort based on the IRGP signature. The signature model was presented as a (Table 2).

What is more, we got the cut-off value of “-0.464” by utilizing the time-dependent ROC curve analysis (Figure 1). According to the cut-off value, the patients were divided into two groups: the high-risk group and low-risk group (Table S1).

(a)

(b)

(a)
(b)

Figure 1 

17-IRGP prognostic signature analysis flowchart (a). Time-dependent ROC analysis of the 17-IRGP prognostic signature for CC (b). Patients were stratified into the high-risk or low-risk groups based on cut-off value “-0.464”, and the AUC of the ROC curve was 0.82.

3.2. Verification of the IRGP Risk Model

To further confirm the prognostic value of the IRGP signature in CC patients, we analyzed the survival of the high-risk and low-risk samples in the discovery cohort and the validation cohort, respectively. Both cohorts suggested that overall survival (OS) in the low-risk group was significantly better than that in the high-risk group (), and the results are shown in Figure 2.

(a)

(b)

(a)
(b)

Figure 2 

Kaplan-Meier curve of OS for CC. According to the IRGP signature, patients were stratified into two groups: the high-risk group (red) and low-risk group (blue). OS of patients in the discovering cohort (a) and validation cohort (b).

For the sake of exploring whether the IRGP signature was independent of other clinical features, we performed Cox regression analysis of univariate and multivariate in two cohorts, and the results are shown in forest plots (Figure 3). Univariate Cox analysis indicated that age, tumor stage, and the IRGP risk score had a significant influence on prognosis in both two datasets (Figures 3(a) and 3(b)).

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 3 

Forest plots of the hazard ratio for assessing the independent prognostic value of the IRGP signature. The univariate (a, b) and multivariate (c, d) Cox regression analysis of age, gender, stage, and risk score in the discovering cohort (a, c) and validation cohort (b, d).

In multivariate Cox regression analysis, the risk score based on IRGP signature was an independent prognostic factor in the discovery cohort (, HR: 3.234, 95% CI: 2.381-4.394) and in the validation cohort (, HR: 1.309, 95% CI: 2.381-4.394). Besides, age and tumor stage were also statistically significant (Figures 3(c) and 3(d)).

3.3. Immune Cell Infiltration Correlated with IRGP Signature

We exploited CIBERSORT to estimate the abundance of 22 kinds of immune cells in each sample of two risk groups in the TCGA dataset. Radar plot showed the distribution of 22 kinds of immune cells in both the high-risk group and low-risk group (Figure 4(a)), among which T cell CD4 memory resting in the low-risk group was significantly higher than that in the high-risk group (); the graphical representation of the result is shown in Figure 4(b).

(a)

(b)

(a)
(b)

Figure 4 

Immune cell infiltration between the two risk groups. (a) The radar plot shows the abundance of 22 immune cells estimated by CIBERSOFT within the two risk groups. (b) The significant distribution of immune cells. The expression of T cell CD4 memory resting was higher in the low-risk group () (, , and ).

3.4. Biological Function Related to the IRGP Signature

To explore the biological roles of differentially expressed genes between the high-risk and low-risk immune groups, we conducted a gene set enrichment analysis on the TCGA dataset. The results showed that genes enriched in seven pathways in the high-risk group, including “cell adhesion molecules cams” and “focal adhesion,” while the expression of genes in the low-risk group was active in three pathways, including “citrate cycle TCA cycle” and “RNA degradation.” The results are shown in Figure 5 and Table S2.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 5 

The biological function associated with the IRGP signature. GSEA shows the CC-related KEGG pathways (). The pathways enriched in the high-risk group (a, b), and the low-risk group (c, d).

4. Discussion

In this study, we used the k-TSP (–top scoring pairs) [19] approach based on the theory of relative expression analysis (RXA) [20] to screen out the immune gene pairs related to colon cancer tissues and constructed a model to predict the prognosis of colon cancer patients according to their characteristics. This method avoids the technical deviation of standardizing gene expression profiles caused by different sequencing platforms [21], and the method can be performed by the web tool (ESurv) [22]. It shows the gene expression value through relative sequencing and pairwise comparison of the same sample, thereby enhancing the robustness of the prediction model.

Herein, the 17 immune gene pairs we screened from the TCGA dataset are composed of 26 immune genes, most of which play the role of antimicrobials, chemokines, and cytokines in the immune process. Among the 26 immune genes, RBP7 has been proved to be highly expressed as an independent biomarker of poor cancer-specific survival in early and advanced CC and is associated with CC progression [23]. The overexpression of CCL4 in CC may induce the infiltration of tumor-associated macrophages, especially the distribution of pretumor macrophages, and there was a positive correlation between plasmatic CCL4 and inflammatory mediators, which have been suggested a poor prognosis [24]. Besides, a previous study indicated that STC2 activates ERK/MEK and PI3K/AKT signaling pathways to promote colorectal tumorigenesis and epithelial-mesenchymal transition progression, and the high expression of STC2 in serum and tumor tissues is related to the low survival rate [25]. In addition, overexpression of miR-4709 promotes the proliferation and invasion of CC by downregulating NR3C2, which is an unfavorable prognostic factor [26].

Meanwhile, our model suggested CC patients with low immune risk had a better prognosis, which was subsequently verified in the validation cohort. The results showed that the content of T cell CD4 memory resting in the high immune-risk group was significantly lower than that in the low immune-risk group, which acts as the protective factor [27, 28].

The results of GSEA indicated that the differentially expressed immune-related genes between the two risk groups were enriched in the pathways: “citrate cycle TCA cycle,” “RNA degradation,” “focal adhesion,” and “cell adhesion molecules (CAMs).” Therein, the “RNA degradation” pathway can regulate the function of lncRNA GAS5 in mammalian cells [29], while GAS5 contributes to not only the susceptibility but also the lymphatic metastasis of colorectal cancer [30]. The other three pathways have been proved to be related to the tumorigenesis and progression of CC [31–33].

Although the present study merged the immune and prognostic features of colon cancer and accurately stratified the prognosis of patients, this model still has drawbacks. First of all, this study is a retrospective analysis of public datasets, which is susceptible to bias, so a large sample prospective study is needed to verify the stability of the signature model; secondly, the model based on gene-level prognostic features to predict the cost of samples is high, and the resistance of clinical promotion is relatively large.

5. Conclusion

In summary, we conducted a comprehensive analysis of the prognostic value of immune-related gene pairs and constructed a model that can provide a risk assessment for the treatment of colon cancer patients, so as to help them benefit from immunotherapy.

Abbreviations

Data Availability

The data used to support the findings of this study is included within the article.

Conflicts of Interest

The authors declare that they have no conflicts.

Authors’ Contributions

Qianshi Zhang and Zhen Feng contributed equally to this work.

Acknowledgments

This work was supported by grants from the Liaoning Province Nature Science Foundation of China (No. 2019-ZD-0917) and the Basic Research Program of Higher Education of Liaoning Province of China (No. LQ2017035).

Supplementary Materials

Supplementary Table 1: Patients’ immune risk stratification. Supplementary Table 2: GSEA. (Supplementary Materials)

References

1. R. L. Siegel, K. D. Miller, A. G. Sauer et al., “Colorectal cancer statistics, 2020,” CA: A Cancer Journal for Clinicians, vol. 70, no. 3, pp. 145–164, 2020. View at: Publisher Site | Google Scholar
2. F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, and A. Jemal, “Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries,” CA: A Cancer Journal for Clinicians, vol. 68, no. 6, pp. 394–424, 2018. View at: Publisher Site | Google Scholar
3. P. Favoriti, G. Carbone, M. Greco, F. Pirozzi, R. E. Pirozzi, and F. Corcione, “Worldwide burden of colorectal cancer: a review,” Updates in Surgery, vol. 68, no. 1, pp. 7–11, 2016. View at: Publisher Site | Google Scholar
4. D. X. Yang, C. P. Gross, P. R. Soulos, and J. B. Yu, “Estimating the magnitude of colorectal cancers prevented during the era of screening: 1976 to 2009,” Cancer, vol. 120, no. 18, pp. 2893–2901, 2014. View at: Publisher Site | Google Scholar
5. T. Hu, Z. Li, C. Y. Gao, and C. H. Cho, “Mechanisms of drug resistance in colon cancer and its therapeutic strategies,” World Journal of Gastroenterology, vol. 22, no. 30, pp. 6876–6889, 2016. View at: Publisher Site | Google Scholar
6. K. Van der Jeught, H.-C. Xu, Y.-J. Li, X.-B. Lu, and G. Ji, “Drug resistance and new therapies in colorectal cancer,” World Journal of Gastroenterology, vol. 24, no. 34, pp. 3834–3848, 2018. View at: Publisher Site | Google Scholar
7. M. Yamadera, E. Shinto, Y. Kajiwara et al., “Differential survival benefits of 5-fluorouracil-based adjuvant chemotherapy for patients with microsatellite-stable stage III colorectal cancer according to the tumor budding status: a retrospective analysis,” Diseases of the Colon and Rectum, vol. 62, no. 11, pp. 1316–1325, 2019, Epub 2019/10/01. View at: Publisher Site | Google Scholar
8. D. R. Almquist, D. H. Ahn, and T. S. Bekaii-Saab, “The role of immune checkpoint inhibitors in colorectal adenocarcinoma,” BioDrugs, vol. 34, no. 3, pp. 349–362, 2020. View at: Publisher Site | Google Scholar
9. D. Basile, S. K. Garattini, M. Bonotto et al., “Immunotherapy for colorectal cancer: where are we heading?” Expert Opinion on Biological Therapy, vol. 17, no. 6, pp. 709–721, 2017. View at: Publisher Site | Google Scholar
10. J. Couzin-Frankel, “Breakthrough of the year 2013. Cancer immunotherapy,” Science., vol. 342, no. 6165, pp. 1432-1433, 2013. View at: Publisher Site | Google Scholar
11. K. M. Bever and D. T. Le, “An expanding role for immunotherapy in colorectal cancer,” Journal of the National Comprehensive Cancer Network, vol. 15, no. 3, pp. 401–410, 2017. View at: Publisher Site | Google Scholar
12. Y. Yang, “Cancer immunotherapy: harnessing the immune system to battle cancer,” Journal of Clinical Investigation, vol. 125, no. 9, pp. 3335–3337, 2015. View at: Publisher Site | Google Scholar
13. L. Marisa, A. de Reyniès, A. Duval et al., “Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value,” PLoS Medicine, vol. 10, no. 5, article e1001453, 2013. View at: Publisher Site | Google Scholar
14. A. Stein and G. Folprecht, “Immunotherapy of colon cancer,” Oncology Research and Treatment, vol. 41, no. 5, pp. 282–285, 2018. View at: Publisher Site | Google Scholar
15. N. J. Llosa, M. Cruise, A. Tam et al., “The vigorous immune microenvironment of microsatellite instable colon cancer is balanced by multiple counter-inhibitory checkpoints,” Cancer Discovery, vol. 5, no. 1, pp. 43–51, 2015. View at: Publisher Site | Google Scholar
16. A. Passardi, M. Canale, M. Valgiusti, and P. Ulivi, “Immune checkpoints as a target for colorectal cancer treatment,” International Journal of Molecular Sciences, vol. 18, no. 6, p. 1324, 2017. View at: Publisher Site | Google Scholar
17. C. Zou, H. Qiu, W. Tang, Y. Wang, B. Lan, and Y. Chen, “CTLA4 tagging polymorphisms and risk of colorectal cancer: a case-control study involving 2,306 subjects,” OncoTargets and Therapy, vol. 11, pp. 4609–4619, 2018. View at: Publisher Site | Google Scholar
18. A. M. Newman, C. L. Liu, M. R. Green et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nature Methods, vol. 12, no. 5, pp. 453–457, 2015. View at: Publisher Site | Google Scholar
19. A. C. Tan, D. Q. Naiman, L. Xu, R. L. Winslow, and D. Geman, “Simple decision rules for classifying human cancers from gene expression profiles,” Bioinformatics, vol. 21, no. 20, pp. 3896–3904, 2005. View at: Publisher Site | Google Scholar
20. J. A. Eddy, J. Sung, D. Geman, and N. D. Price, “Relative expression analysis for molecular cancer diagnosis and prognosis,” Technology in Cancer Research & Treatment, vol. 9, no. 2, pp. 149–159, 2010. View at: Publisher Site | Google Scholar
21. A. S. Devonshire, R. Elaswarapu, and C. A. Foy, “Evaluation of external RNA controls for the standardisation of gene expression biomarker measurements,” BMC Genomics, vol. 11, no. 1, p. 662, 2010. View at: Publisher Site | Google Scholar
22. K. Pak, S. O. Oh, T. S. Goh et al., “A user-friendly, web-based integrative tool (ESurv) for survival analysis: development and validation study,” Journal of Medical Internet Research, vol. 22, no. 5, article e16084, 2020. View at: Publisher Site | Google Scholar
23. M. Elmasry, L. Brandl, J. Engel, A. Jung, T. Kirchner, and D. Horst, “RBP7 is a clinically prognostic biomarker and linked to tumor invasion and EMT in colon cancer,” Journal of Cancer, vol. 10, no. 20, pp. 4883–4891, 2019. View at: Publisher Site | Google Scholar
24. M. De la Fuente López, G. Landskron, D. Parada et al., “The relationship between chemokines CCL2, CCL3, and CCL4 with the tumor microenvironment and tumor-associated macrophage markers in colorectal cancer,” Tumor Biology, vol. 40, no. 11, 2018. View at: Publisher Site | Google Scholar
25. B. Chen, X. Zeng, Y. He et al., “STC2 promotes the epithelial-mesenchymal transition of colorectal cancer cells through AKT-ERK signaling pathways,” Oncotarget, vol. 7, no. 44, pp. 71400–71416, 2016. View at: Publisher Site | Google Scholar
26. M. Yu, H. L. Yu, Q. H. Li, L. Zhang, and Y. X. Chen, “miR-4709 overexpression facilitates cancer proliferation and invasion via downregulating NR3C2 and is an unfavorable prognosis factor in colon adenocarcinoma,” Journal of Biochemical and Molecular Toxicology, vol. 33, no. 12, article e22411, 2019. View at: Publisher Site | Google Scholar
27. P. Ge, W. Wang, L. Li et al., “Profiles of immune cell infiltration and immune-related genes in the tumor microenvironment of colorectal cancer,” Biomedicine & Pharmacotherapy, vol. 118, p. 109228, 2019. View at: Publisher Site | Google Scholar
28. X. Zhang, F. Quan, J. Xu, Y. Xiao, X. Li, and Y. Li, “Combination of multiple tumor-infiltrating immune cells predicts clinical outcome in colon cancer,” Clinical Immunology, vol. 215, p. 108412, 2020. View at: Publisher Site | Google Scholar
29. H. Tani, M. Torimura, and N. Akimitsu, “The RNA degradation pathway regulates the function of GAS5 a non-coding RNA in mammalian cells,” PLoS One, vol. 8, no. 1, p. e55684, 2013. View at: Publisher Site | Google Scholar
30. Y. Zheng, D. Song, K. Xiao et al., “LncRNA GAS5 contributes to lymphatic metastasis in colorectal cancer,” Oncotarget, vol. 7, no. 50, pp. 83727–83734, 2016. View at: Publisher Site | Google Scholar
31. M. Kitazawa, T. Hatta, Y. Sasaki et al., “Promotion of the Warburg effect is associated with poor benefit from adjuvant chemotherapy in colorectal cancer,” Cancer Science, vol. 111, no. 2, pp. 658–666, 2020. View at: Publisher Site | Google Scholar
32. B. T. Nguyen, J. C. Pyun, S. G. Lee, and M. J. Kang, “Identification of new binding proteins of focal adhesion kinase using immunoprecipitation and mass spectrometry,” Scientific Reports, vol. 9, no. 1, article 12908, 2019. View at: Publisher Site | Google Scholar
33. X. Bessa, J. I. Elizalde, F. Mitjans et al., “Leukocyte recruitment in colon cancer: Role of cell adhesion molecules, nitric oxide, and transforming growth factor β1,” Gastroenterology, vol. 122, no. 4, pp. 1122–1132, 2002. View at: Publisher Site | Google Scholar

Copyright

Copyright © 2021 Qianshi Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.