ORIGINAL RESEARCH article

Front. Genet., 04 March 2021

Sec. Cancer Genetics

Volume 12 - 2021 | https://doi.org/10.3389/fgene.2021.619611

Prognostic Risk Model of Immune-Related Genes in Colorectal Cancer

  • 1. Department of Colorectal Surgery and Oncology, Key Laboratory of Cancer Prevention and Intervention, Ministry of Education, The Second Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China

  • 2. Zhejiang University Cancer Center, Zhejiang University, Hangzhou, China

Abstract

Purpose:

We focused on immune-related genes (IRGs) derived from transcriptomic studies, which had the potential to stratify patients’ prognosis and to establish a risk assessment model in colorectal cancer.

Summary:

This article examined our understanding of the molecular pathways associated with intratumoral immune response, which represented a critical step for the implementation of stratification strategies toward the development of personalized immunotherapy of colorectal cancer. More and more evidence shows that IRGs play an important role in tumors. We have used data analysis to screen and identify immune-related molecular biomarkers of colon cancer. We selected 18 immune-related prognostic genes and established models to assess prognostic risks of patients, which can provide recommendations for clinical treatment and follow-up. Colorectal cancer (CRC) is a leading cause of cancer-related death in human. Several studies have investigated whether IRGs and tumor immune microenvironment (TIME) could be indicators of CRC prognoses. This study aimed to develop an improved prognostic signature for CRC based on IRGs to predict overall survival (OS) and provide new therapeutic targets for CRC treatment. Based on the screened IRGs, the Cox regression model was used to build a prediction model based on 18-IRG signature. Cox regression analysis revealed that the 18-IRG signature was an independent prognostic factor for OS in CRC patients. Then, we used the TIMER online database to explore the relationship between the risk scoring model and the infiltration of immune cells, and the results showed that the risk model can reflect the state of TIME to a certain extent. In short, an 18-IRG prognostic signature for predicting CRC patients’ survival was firmly established.

Introduction

Colorectal cancer (CRC) ranks among the top causes of cancer-related deaths worldwide that endangers human health. The GLOBOCAN data in 2018 released by the International Cancer Research Agency showed that each year there were approximately 1.85 million new CRCs and more than 880,000 deaths worldwide. The morbidity and mortality of CRC rank third and second, respectively, in malignant tumors, in which the morbidity accounts for approximately 10% of the total cancer incidence, and the mortality accounts for 9% of the total deaths due to cancer (). It was predicted that the number of cases will increase by more than 60% in 2030, with 2.2 million new cases and 1.1 million deaths (). Surgical resection is the main treatment option for CRC patients. With the application and popularity of colonoscopy, early treatment work has been improved. The clinical outcomes of CRC patients in many countries have improved significantly over the past few decades (). Despite the complete surgical resection, many CRC patients eventually relapsed and developed metastatic disease (). In clinical practice, a more effective prognostic evaluation system is urgently needed to provide personalized medicine for CRC patients and improve patient outcomes.

It is noteworthy that after Fearon and Vogelstein proposed the model of CRC genetic basis, researchers have begun to understand the heterogeneity of CRC (). Patients with different genetic backgrounds had different outcomes after receiving the same treatment (). Some researchers believed that it was attributed to immunity-related factors (). As we knew that the immune system plays an important role in the development of a variety of cancers, including CRC (). A recent study found that immunological data (such as type, density, and location of immune cells in tumor samples) can predict patient survival better than the current histopathological characteristics used for CRC patients (). Immune cells are important parts of the tumor microenvironment and affect the development and metastasis of CRC (). Tumor-infiltrating macrophages and dendritic cells in CRC are related to local regulatory T cells and systemic T-cell responses to tumor-associated antigens and have an impact on patients’ survival (). In addition, studies have shown that immune-related genes (IRGs) in colon cancer are closely related to the occurrence and development of colon cancer. However, there is currently no prognostic model based on IRGs to predict the overall prognosis of CRC patients and systematically assess the immune environment of CRC (). Therefore, constructing an immune-based prognostic model that can effectively predict the prognosis of CRC has a very important clinical application prospect.

In this study, we screened differentially expressed IRGs that are closely related to CRC through bioinformatics analysis of The Cancer Genome Atlas (TCGA). Next, the IRGs that were significantly associated with prognosis were further screened. Differentially expressed tumor-associated transcription factors (TFs) were searched, and a correlation network was constructed to reveal the relationship between TFs regulating immune genes. Then, immune-related prognostic models were constructed by integrating IRGs of CRC. Besides, we verified that the risk model can be used as an effective independent prognostic indicator.

Materials and Methods

Patient Data Collection

Colorectal cancer patients (adenocarcinomas) with gene expression profiles and clinical information were obtained from TCGA data portal1. Processed RNA-Seq FPKM data of 398 CRC and 39 adjacent normal tissues were downloaded for further analyses.

IRGs and Cancer-Related Transcription Factors

The comprehensive list of IRGs was downloaded from the Immunology Database and Analysis Portal (ImmPort) database2, which shares immunology data and provides a list of IRGs for cancer researchers (). The IRGs that actively participated in the immune process were identified. To investigate the regulatory mechanism of IRGs, we extracted cancer-related transcription factors (CRTFs) for subsequent research. The CRTF data were downloaded from the Cistrome Cancer database3, which is a useful database for biomedical and genetic research and includes 318 CRTFs ().

Differential Gene Expression Analysis

To select the IRGs and TFs that contributed to the development and progression of CRC, differentially expressed genes (DEGs) between tumor samples and normal samples were screened using the limma R package. Differential expression analysis was conducted, with an adjusted false discovery rate < 0.05 and | log2(fold change)| > 1 as the thresholds. Differentially expressed IRGs were identified as overlaps between the IRG list and the DEG list. Differentially expressed TFs were identified as overlaps between the TFs list and the DEG list. Heatmaps were generated using the “pheatmap” R package, and volcano plots were also displayed using the “ggplot2” R package.

CRTF-IRG Regulatory Network

In order to evaluate how differentially expressed CRTFs regulate prognosis-related IRGs, we studied the correlation between them. The core method is the Pearson test. The critical standard is set to a correlation coefficient > 0.4, P < 0.001. This step is performed using the Cor. test function in R, and the correlation coefficient and P value are calculated by Cor. test. To make the situation clearer, Cytoscape was used to build a visual regulatory network.

PPI Network Construction and Module Analysis

The PPI network was predicted using an online database search tool STRING4 (). Analyzing functional interactions between proteins may provide insights into the mechanisms of CRC development and progression. In this study, prognostic-related PPI networks of IRGs and CRTFs were constructed using the STRING database, and interactions with composite scores > 0.4 were considered statistically significant. Kyoto Encyclopedia of Genes and Genomes signaling pathways and biological functions of genes were analyzed using functional clustering carried by STRING.

Gene Set Enrichment Analysis

Gene Set Enrichment Analysis (GSEA)5 was used to analyze the GO term of the genes that make up the signature.

Construction of the Immune-Related Signature for CRC

To control the quality of the data, after excluding patients who lacked survival information or survived for less than 90 days, 334 samples were subsequently analyzed. Transcriptomic analysis of RNA measured by FPKM values was performed using log2-based conversion. Based on the differentially expressed IRGs, Kaplan–Meier analysis was first performed to screen prognostic immune genes. Then, immune-related prognostic signature (IPS) was constructed by multivariate Cox regression to calculate the risk score for each patient. Risk scores were acquired based on expressions of genes multiplied by a linear combination of regression coefficients obtained from the multivariate Cox regression analysis. P < 0.01 was regarded as significant.

Survival Analysis

According to the optimal cutoff value obtained by the “survminer” R package, CRC patients were classified into low risk and high risk according to their risk scores. To investigate the prognostic value of the prognostic model in CRC patients, univariate Cox analysis was implemented by the “survival” R package and “survminer” R package. A time-dependent receiver operating characteristic (ROC) curve was plotted to assess sensitivity and specificity using the “timeROC” R software package (). The area under the curve was calculated from the ROC curve.

Association Analysis Between 18-IRGs and Clinical Parameters

Association analysis of clinical characteristics of 18 key prognostic IRGs in the model was performed using the t test. To transform the data types into binary variables, 398 CRC patients were grouped according to different clinical characteristics. In terms of age, 65 years old was chosen as the cutoff point. The stage was divided into stages I and II and stages III and IV. The T stage was divided into T1–2 and T3–4. M stage was divided into M0 or M1. N stages N0 and N1–2.

TIMER Database Analysis of the Correlation Between Immune-Related Markers and Immune Cell Infiltration

TIMER database6 is a comprehensive resource for systematical analysis of immune infiltrates across different cancer types (). The abundance of six immune infiltrates was estimated by the TIMER algorithm (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells). We used the TIMER database to analyze the correlation between the prognostic model of CRC patients and six tumor-infiltrating immune cells.

Statistical Analysis

Overall survival (OS) was defined as the main outcome. Univariate cox regression analysis and multivariate cox regression analysis were performed to evaluate the prognostic effect of the immune signature and various clinicopathological features including age, clinical stage, grade, and TNM stage. Statistical analyses were performed using R software (version 3.5.1). The heatmap was generated using the “pheatmap” R package. Unless otherwise specified, a two-sided P < 0.05 was considered statistically significant.

Results

Differentially Expressed IRGs and CRTFs in CRC

Compared with normal tissues, there were 5,938 DEGs in CRC tissues, of which 3,936 were up-regulated and 2,002 were down-regulated in these samples. The difference between tumor tissue and normal tissue can be seen through the heatmap and the volcano map (Figures 1A,B). Compared with normal tissues, a total of 484 IRGs (173 up-regulated and 311 down-regulated) and 71 CRTFs (46 up-regulated and 25 down-regulated) were differentially expressed in CRC tissues. The heatmaps showed that CRC samples can be distinguished from normal samples based on the differentially expressed IRGs and CRTFs (Figures 1C,D). The volcano plots showed the distribution of differentially expressed IRGs and CRTF between CRC samples and normal controls (Figures 1E,F).

FIGURE 1

Screening of IRGs Related to Significant Prognosis in CRC

To determine the differentially expressed IRGs with prognostic characteristics, the relationship between the expression of 484 IRGs in 398 CRC samples and prognosis were evaluated by univariate Cox analysis. A total of 30 IRGs with prognostic characteristics were found, as shown in Table 1. Figure 2 is a forest plot showing the prognostic IRGs, P values, and hazard ratios. Among the 18 prognostic-related IRGs, CD1B, CXCL3, F2RL1, and IGHG4 are low-risk genes. The higher expression of these genes indicated better prognosis of patients. The other 14 IRGs are high-risk genes, and when their expression increases, the patient’s risk increases. NGF is the gene with the highest risk factor.

TABLE 1

Gene symbolHRHR.95LHR.95HP-valueGene symbolHR (95% CI)P-value
CD1B0.0570.0070.4690.007658134CD1B0.057 (0.007–0.469)0.008
SLC10A21.8331.1742.8600.007661263SLC10A21.833 (1.174–2.860)0.008
CXCL30.9760.9600.9940.007314764CXCL30.976 (0.960–0.994)0.007
NOX41.6441.1612.3280.00512998NOX41.644 (1.161–2.328)0.005
FABP41.0141.0071.0204.75E-05FABP41.014 (1.007–1.020)4.75E-05
ADIPOQ1.1021.0471.1600.00022247ADIPOQ1.102 (1.047–1.160)2.22E-04
FGF21.3401.0761.6700.009000359FGF21.340 (1.076–1.670)0.009
F2RL10.9670.9440.9910.006461334F2RL10.967 (0.944–0.991)0.006
CCL191.0301.0081.0510.005852187CCL191.030 (1.008–1.051)0.006
PLCG21.6741.1922.3510.002941416PLCG21.674 (1.192–2.351)0.003
IGHG11.0011.0001.0010.000824658IGHG11.001 (1.000–1.001)0.001
IGHG41.0001.0001.0010.008255942IGHG41.000 (1.000–1.001)0.008
IGHV4-311.0081.0021.0140.007667535IGHV4-311.008 (1.002–1.014)0.008
IGHV5-511.0021.0011.0030.002062246IGHV5-511.002 (1.001–1.003)0.002
IGKV1-331.0301.0111.0500.001978816IGKV1-331.030 (1.011–1.050)0.002
IGKV1-81.0441.0151.0730.002316247IGKV1-81.044 (1.015–1.073)0.002
IGKV2D-401.0161.0051.0260.003078753IGKV2D-401.016 (1.005–1.026)0.003
IGLV6-571.0021.0011.0030.003346154IGLV6-571.002 (1.001–1.003)0.003
SEMA3G1.2941.1231.4910.000355359SEMA3G1.294 (1.123–1.491)3.55E-04
INHBA1.0531.0221.0850.000678328INHBA1.053 (1.022–1.085)0.001
NGF3.6152.0386.4131.12E-05NGF3.615 (2.038–6.413)1.12E-05
RETNLB1.0041.0011.0060.002823016RETNLB1.004 (1.001–1.006)0.003
STC11.0781.0211.1390.007200166STC11.078 (1.021–1.139)0.007
UCN1.3831.1181.7110.002826168UCN1.383 (1.118–1.711)0.003
VIP1.0581.0211.0960.001913177VIP1.058 (1.021–1.096)0.002
NGFR1.2001.0921.3200.000160985NGFR1.200 (1.092–1.320)1.61E-04
NPR11.5011.1331.9880.004630474NPR11.501 (1.133–1.988)0.005
OXTR1.4261.1771.7280.000284786OXTR1.426 (1.177–1.728)2.85E-04
PTH1R1.6281.2202.1740.000936214PTH1R1.628 (1.220–2.174)0.001
TRDC1.1491.0361.2740.008408965TRDC1.149 (1.036–1.274)0.008

General characteristics of prognostic immune-related genes.

FIGURE 2

The Mechanism of Prognosis-Related IRGs and CRTF-IRG Regulatory Network

We explored the potential regulatory mechanisms of 18 prognostic-related IRGs, which may reflect the regulatory mechanisms of these gene sets. We selected 30 prognostic-related IRGs and 71 differential CTRFs for correlation analysis to explore the regulatory mechanism of prognostic-related IRGs. The Cor. test function is used to test the correlation between each CRTF and each IRG. The core method is Pearson test. The correlation coefficient filter is 0.4, and the P value filter is 0.001. The regulatory relationship between these CRTFs and IRGs is revealed in the regulatory network (Figure 3A). As shown in Figure 3A, NR3C1, MYH11, RUNX1, MAF, CCB7, LMO2, FOXP3, and EPAS1 regulate most of the IRGs related to prognosis and dominate the regulation network. This transcriptional regulatory network reveals the regulatory relationship between these IRGs and CRTFs. Table 2 shows the correlation between IRGs and CRTFs after screening. The PPI network of IRGs and CRTFs was constructed, and the most significant module was obtained. The functional analyses of genes involved in this module were analyzed. Enrichment analysis shows that the genes in this module are mainly involved in cell proliferation and metabolic processes (Figure 3B).

FIGURE 3

TABLE 2

TFImmune GeneCorP valueRegulation
BHLHE40CLCF10.4168393141.80E-15Positive
BHLHE40INHBA0.4015073852.28E-14Positive
CBX7CXCL120.4865617662.97E-21Positive
CBX7PTGDS0.5032001597.70E-23Positive
CBX7COLEC120.4240782955.19E-16Positive
CBX7A2M0.6162427872.60E-36Positive
CBX7CCL190.4365622945.65E-17Positive
CBX7SEMA3G0.5427597565.56E-27Positive
CBX7SLIT20.5558544371.78E-28Positive
CBX7TNFSF120.5632640842.36E-29Positive
CBX7NGFR0.4886033111.92E-21Positive
CBX7NPR10.5312604511.01E-25Positive
CBX7S1PR10.5925046314.92E-33Positive
CDK2BIRC50.4314052351.43E-16Positive
CENPABIRC50.6421776643.16E-40Positive
E2F3S100P−0.443391071.61E-17Negative
EPAS1CXCL120.4537007152.31E-18Positive
EPAS1PTGDS0.4480129056.81E-18Positive
EPAS1A2M0.556190751.62E-28Positive
EPAS1CCL190.4131499743.35E-15Positive
EPAS1PLCG20.4463883369.24E-18Positive
EPAS1SEMA3G0.4841293214.99E-21Positive
EPAS1NPR10.4245772774.75E-16Positive
EPAS1S1PR10.5622560833.12E-29Positive
EZH2BIRC50.4045145181.40E-14Positive
FOSL1CLCF10.5399651721.14E-26Positive
FOXP3CD1B0.5125154439.11E-24Positive
FOXP3PTGDS0.4275768452.81E-16Positive
FOXP3A2M0.5456659372.62E-27Positive
FOXP3TLR70.5059431614.13E-23Positive
FOXP3PLCG20.4448091831.24E-17Positive
FOXP3IGHG10.4453391741.12E-17Positive
FOXP3CMKLR10.650445521.48E-41Positive
FOXP3TNFSF120.4620214034.58E-19Positive
FOXP3S1PR10.4915695631.01E-21Positive
H2AFXBIRC50.4592903017.83E-19Positive
KAT2BNR3C20.4537927522.27E-18Positive
KLF4CCL280.4280002082.61E-16Positive
KLF4GUCA2A0.4379377354.40E-17Positive
KLF4NR3C20.5656588671.22E-29Positive
LMO2CXCL120.5745742379.84E-31Positive
LMO2PTGDS0.5391333081.40E-26Positive
LMO2COLEC120.4104595255.26E-15Positive
LMO2A2M0.5899181061.08E-32Positive
LMO2CCL190.4529420432.67E-18Positive
LMO2CD79B0.435469196.88E-17Positive
LMO2PLCG20.4963048113.59E-22Positive
LMO2SEMA3G0.5750298168.63E-31Positive
LMO2SLIT20.4611475645.44E-19Positive
LMO2CMKLR10.4650547832.51E-19Positive
LMO2NGF0.4086958767.04E-15Positive
LMO2TNFSF120.4975801812.70E-22Positive
LMO2NGFR0.4563161871.40E-18Positive
LMO2NPR10.5468690581.92E-27Positive
LMO2S1PR10.6353951463.63E-39Positive
MAFCXCL120.6604488843.22E-43Positive
MAFPTGDS0.5247812964.95E-25Positive
MAFCOLEC120.7008042511.23E-50Positive
MAFA2M0.720798678.62E-55Positive
MAFNOX40.6092113622.60E-35Positive
MAFTLR70.6584408767.02E-43Positive
MAFPLCG20.5072154463.09E-23Positive
MAFIGHG10.4370030345.21E-17Positive
MAFSEMA3G0.5816102921.28E-31Positive
MAFSLIT20.6212766314.82E-37Positive
MAFCMKLR10.6836283252.48E-47Positive
MAFINHBA0.597944249.22E-34Positive
MAFNGF0.4301768791.78E-16Positive
MAFSTC10.4386396683.87E-17Positive
MAFTNFSF120.6322654461.10E-38Positive
MAFNGFR0.4177314111.55E-15Positive
MAFNPR10.5908357078.17E-33Positive
MAFS1PR10.7147783521.68E-53Positive
MYH11CXCL120.409253576.42E-15Positive
MYH11PTGDS0.4223867926.95E-16Positive
MYH11A2M0.6137263695.96E-36Positive
MYH11SEMA3G0.4276127392.79E-16Positive
MYH11SLIT20.5558333481.79E-28Positive
MYH11TNFSF120.4615761075.00E-19Positive
MYH11VIP0.4704637488.48E-20Positive
MYH11NGFR0.4561879731.43E-18Positive
MYH11NPR10.5204535081.40E-24Positive
MYH11S1PR10.5900453061.04E-32Positive
NCAPGBIRC50.5548397182.33E-28Positive
NR3C1CXCL120.5890989231.38E-32Positive
NR3C1PTGDS0.4495037035.14E-18Positive
NR3C1COLEC120.625468251.16E-37Positive
NR3C1A2M0.6753231038.15E-46Positive
NR3C1NOX40.5023867999.25E-23Positive
NR3C1TLR70.5490007061.10E-27Positive
NR3C1FGF20.4567142871.29E-18Positive
NR3C1PLCG20.4652573172.41E-19Positive
NR3C1SEMA3G0.5120801431.01E-23Positive
NR3C1SLIT20.5842750415.83E-32Positive
NR3C1CMKLR10.5478170161.50E-27Positive
NR3C1INHBA0.5596006096.44E-29Positive
NR3C1TNFSF120.4823806737.22E-21Positive
NR3C1VIP0.4207462469.23E-16Positive
NR3C1IL1RAP0.4097895765.88E-15Positive
NR3C1NPR10.468519971.26E-19Positive
NR3C1S1PR10.6461790397.28E-41Positive
PBX1A2M0.4050348891.28E-14Positive
RUNX1CXCL120.4382092684.19E-17Positive
RUNX1COLEC120.4075396718.52E-15Positive
RUNX1A2M0.4333017181.02E-16Positive
RUNX1SEMA3G0.4165834171.88E-15Positive
RUNX1INHBA0.4358417246.44E-17Positive
RUNX1S1PR10.402023122.10E-14Positive
SNAPC4JAG20.4022111182.03E-14Positive
SOX4S100P−0.415710928  2.18E-15Negative
SPDEFS100P0.4250639764.37E-16Positive
SPDEFRETNLB0.4958661863.95E-22Positive
SPIBSCG20.4198921771.07E-15Positive
TFAP2CIGHV3-640.5077063922.76E-23Positive
TFAP2CPTH1R0.4592182827.94E-19Positive

Correlation between prognostic IRGs and CRTFs.

Hub Gene Selection and Analysis in CRC

Using Cytohubba in Cytoscape, we filtered 33 hub genes that were identified by filtering according to the criterion of degrees > 10 criteria (each node had more than 10 interactions), and the 10 most central genes in the immune gene regulatory network according to node degree were MAF, A2M, CBX7, MYH11, EPAS1, CXCL12, LMO2, S1PR1, FOXP3, and NR3C1 (Figures 4A,B). Gene Ontology (GO) enrichment analysis of genes in the immune gene regulatory network related to prognosis was conducted to explore which signaling pathways were activated. The analysis of the biological processes (BPs) of the central genes using BiNGO in Cytoscape is shown in Figures 4C,D. GO analysis showed that the changes in the BPs of these genes were significantly enriched in the immune process, cell proliferation, immune organ development, and hemopoiesis. Changes in molecular function were mainly focused on TF activity, cytokine activation, and molecular binding. Afterward, the functional enrichment analysis of the key genes of the IRG set was performed by GSEA. Figure 4E shows that the changes in the BP of these genes are significantly enriched in the immune process, cell proliferation, immune organ development, and hematopoiesis. The results of Kaplan–Meier analysis of these hub genes are in the Supplementary Figure 1.

FIGURE 4

Construction of the Immune-Related Signature for CRC

Multivariate Cox analysis was performed on 30 prognostic IRGs, and 18 genes were finally selected to establish a prognostic model (Table 3). The risk score is based on the gene expression level multiplied by its corresponding regression coefficient. The regression coefficient was calculated by multivariate Cox regression. The risk score is related to not only the expression level of these genes but also the correlation coefficients. The risk score of each patient is the sum of all the 18 risk prognostic genes in Table 3 multiplied by the corresponding risk factors. The 398 CRC samples were then divided into high-risk groups (n = 199) and low-risk groups (n = 199) based on the median risk score (Figure 5A). Survival overview and gene expression heatmaps are presented in Figures 5B,C. Survival analysis showed that the OS of patients in the high-risk group was significantly lower than that in the low-risk group (P < 0.0001; Figure 5D). The 5-year survival rate of the high-risk group was 51.1%, and the 5-year survival rate of the low-risk group was 81.4%. The areas under the ROC curves at 1, 3, and 5 years of OS are 0.811, 0.711, and 0.734, respectively, which indicated that the prognostic model showed good sensitivity and specificity (Figure 5E). In addition, as shown in Supplementary Figure 2, the model after excluding genes with P ≥ 0.05 has advantages in the short-term prognosis (1 year), but the model is not effective in predicting the long-term prognosis.

TABLE 3

Gene symbolCoefHRHR.95LHR.95HP valueGene symbolCoefHR (95% CI)P-value
CD1B−4.725570.0088660.0006070.129560.000554CD1B−4.7260.009 (0.001–0.130)0.001
SLC10A20.8443782.3265291.3931753.8851820.00125SLC10A2  0.8442.327 (1.393–3.885)0.001
CXCL3−0.018820.9813520.9627691.0002940.053627CXCL3−0.0190.981 (0.963–1.000)0.054
NOX4−1.253480.285510.0770171.0584130.060787NOX4−1.2530.286 (0.077–1.058)0.061
FABP40.0569391.0585911.0171121.1017620.00524FABP4  0.0571.059 (1.017–1.102)0.005
ADIPOQ−0.249290.7793530.6157160.9864790.038156ADIPOQ−0.2490.779 (0.616–0.986)0.038
F2RL1−0.026710.973640.9474421.0005610.054904F2RL1−0.0270.974 (0.947–1.001)0.055
PLCG20.4993771.6476951.0490152.5880450.030183PLCG2  0.4991.648 (1.049–2.588)0.030
IGKV1-330.0511571.0524881.0130341.0934790.008684IGKV1-33  0.0511.052 (1.013–1.093)0.009
IGLV6-570.0029351.0029391.0013971.0044840.000186IGLV6-57  0.0031.003 (1.001–1.004)<0.001  
INHBA0.1393991.1495821.0326991.2796950.010831INHBA  0.1391.150 (1.033–1.280)0.011
NGF0.9438962.5699740.9198857.1799920.071757NGF  0.9442.570 (0.920–7.180)0.072
RETNLB0.0041241.0041321.0013741.0068970.003296RETNLB  0.0041.004 (1.001–1.007)0.003
UCN0.4680881.5969381.2511162.0383480.00017UCN  0.4681.597 (1.251–2.038)<0.001  
VIP0.0665151.0687771.0044821.1371870.035621VIP  0.0671.069 (1.004–1.137)0.036
NGFR−0.436370.6463760.4126951.0123760.05662NGFR−0.4360.646 (0.413–1.012)0.057
OXTR−0.304430.7375410.5246611.0367960.079773OXTR−0.3040.738 (0.525–1.037)0.080
TRDC0.2670541.3061111.149421.4841634.21E-05TRDC  0.2671.306 (1.149–1.484)<0.001  

Eighteen genes that constitute the immune-related prognostic model and the corresponding risk factors Riskscore = CD1B*(-4.726) + SLC10A2*(0.844) + CXCL3*(-0.019) + NOX4*(-1.253) + FABP4*(0.057) + ADIPOQ*(-0.249) + F2RL1*(-0.027) + PLCG2*(0.499) + IGKV1 - 33*(0.051) + IGLV6 - 57*(0.003) + INHBA*(0.139) + NGF*(0.944) + RETNLB*(0.004) + UCN*(0.468) + VIP*(0.067) + NGFR*(-0.436) + OXTR*(-0.304) + TRDC*(0.267).

FIGURE 5

Immune-Related Prognostic Signature Was an Independent Predictive Marker of OS for CRC Patients

Three hundred ninety-eight CRC patients with clinical information of age, gender, pathological stage, TNM stage, and risk score were selected for further analysis. Univariate and multivariate Cox regression analyses were performed to assess the independent predictive power of immune-related prognostic markers. Univariate analysis showed that pathological stage (P < 0.001), TNM stage (P < 0.001), and immune-related prognostic risk score (P < 0.001) were significantly correlated with OS (Table 4 and Figure 6A). After multivariate analysis, the immune-related prognostic risk score was the only independent prognostic factor related to OS (P < 0.005; Table 5 and Figure 6B).

TABLE 4

VariableHRHR.95LHR.95HP valueVariableHRP-value
Age1.7360.8963.365  0.102Age1.736 (0.896–3.365)  0.102
Gender1.1780.652.137  0.589Gender1.178 (0.650–2.137)  0.589
Stage2.9082.0394.148<0.001Stage2.908 (2.039–4.148)<0.001
T4.2792.3347.844<0.001T4.279 (2.334–7.844)<0.001
M6.6083.61312.087<0.001M6.608 (3.613–12.087)<0.001
N2.3441.6623.305<0.001N2.344 (1.662–3.305)<0.001
riskScore5.1682.30511.586<0.001riskScore5.168 (2.305–11.586)<0.001

Univariate analyses of overall survival in CRC patients of TCGA.

TABLE 5

VariablesHRHR.95LHR.95HP valueVariablesHR (95% CI)P-value
Age2.3677211.1571184.8448840.018303Age2.368 (1.157–4.845)0.018
Gender1.0942560.5959332.0092810.771426Gender1.094 (0.596–2.009)0.771
Stage1.7618490.6136795.0582040.292554Stage1.762 (0.614–5.058)0.293
T1.6169130.7902573.3083020.188337T1.617 (0.790–3.308)0.188
M1.844770.4521777.5262060.393327M1.845 (0.452–7.526)0.393
N1.1147030.6111022.0333140.723281N1.115 (0.611–2.033)0.723
riskScore3.4729381.4970538.0566940.003734riskScore3.473 (1.497–8.057)0.004

Multivariate analyses of overall survival in CRC patients of TCGA.

FIGURE 6

Association Between 18 IRGs, Clinical Parameters, and Prognostic Risk Scores

We analyzed the association between the expression of 18 key prognostic related IRGs in the patient’s tumor tissue and the patient’s clinical characteristics. The association between NGF, TRDC, CXCL3, CD1B, VIP, F2RL1, FABP4, OXTR, UCN, NOX4, ADIPOQ, and clinical characteristics was found (Table 6 and Figure 7). NGF is negatively correlated with age, and NGF expression is generally higher in advanced patients. Patients with higher VIP expression generally have higher T and N stages. On the other hand, TRDC, CXCL3, and FRL1 are highly expressed in patients in the early stage and patients with N0 stage.

TABLE 6

Gene symbolAgeGenderStageTMN
CD1B−0.983 (0.326)      1.524 (A200.129)2.047 (0.042)  1.837 (0.069)       17.361 (5.956e−04)1.883 (0.061)
SLC10A2     0.57 (0.569)1.211 (0.228)0.673 (0.501)−0.577 (0.565)2.255 (0.521)  0.63 (0.529)
CXCL3−0.593 (0.554)−0.845 (0.399)            3.828 (1.582e−04)  0.582 (0.562)8.458 (0.037)          3.696 (2.634e−04)
NOX4   1.983 (0.048)−0.889 (0.375)  −0.739 (0.460) −0.221 (0.826)1.356 (0.716)−1.146 (0.253) 
FABP4   1.665 (0.098)0.985 (0.326)  −0.97 (0.333)−2.659 (0.008)3.322 (0.345)−1.062 (0.290) 
ADIPOQ     1.49 (0.138)0.232 (0.817)−0.467 (0.641)−2.356 (0.019)1.274 (0.735)−0.578 (0.564) 
F2RL1−1.259 (0.209)0.132 (0.895)  2.675 (0.008)  0.839 (0.404)1.936 (0.586)  2.752 (0.006)
PLCG2−0.344 (0.731)0.971 (0.333)−1.347 (0.179)−0.604 (0.547)2.094 (0.553)−1.693 (0.092)
IGKV1-33−0.892 (0.373)−1.171 (0.243)    1.102 (0.272)−0.602 (0.548)5.168 (0.160)  1.094 (0.275)
IGLV6-57  −0.47 (0.639)−0.851 (0.396)  −0.451 (0.653)  0.023 (0.982)3.862 (0.277)−0.532 (0.596)
INHBA  1.674 (0.095)−0.644 (0.520)  −0.861 (0.390)−0.854 (0.396)1.453 (0.693)−1.192 (0.234)
NGF  1.982 (0.049)   0.97 (0.333)−2.586 (0.010)−2.025 (0.045)3.632 (0.304)−2.864 (0.005)
RETNLB  0.556 (0.579)   0.35 (0.727)  1.355 (0.176)  1.381 (0.172)  1.72 (0.633)  1.273 (0.204)
UCN−2.129 (0.034) 1.217 (0.225)−1.575 (0.117)  −0.26 (0.795)0.735 (0.865)−1.428 (0.155)
VIP  0.486 (0.627)−0.139 (0.889)−1.763 (0.080)−2.259 (0.025)0.968 (0.809)−2.041 (0.043)
NGFR  1.548 (0.124)  1.651 (0.100)−1.511 (0.133)−1.652 (0.101)4.464 (0.216)−1.626 (0.106)
OXTR  0.297 (0.767)−1.763 (0.080)−1.243 (0.215)−1.985 (0.048)1.805 (0.614)−1.272 (0.205)
TRDC  0.144 (0.885)−0.141 (0.888)  2.772 (0.006)   1.011 (0.316)14.881 (0.002)  2.671 (0.008)
riskScore−1.146 (0.253)−0.901 (0.369)−1.268 (0.207)−1.394 (0.165)          17.773 (4.899e−04)−1.274 (0.205)

Eighteen genes in the risk score model and clinical characteristics correlation analysis.

The numbers in the table represent the t value of t test between each gene and clinical features; the numbers in parentheses represent P value.

FIGURE 7

TIMER Database Analysis

The relationships between the risk score model and immune cell infiltration were studied. The characterization of immune infiltration is very important for exploring the state of the immune microenvironment and studying the interaction between tumors and immunity. We applied the TIMER tool to identify potential relationships between IPS and infiltrating immune cells, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. As shown in Figure 8, the proportions of tumor-infiltrating CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells were closely related to our prognostic risk score (p < 0.05).

FIGURE 8

Discussion

In recent years, the genetic characteristics of mRNA in cancer patients have attracted people’s attention, and studies have revealed its great potential in the prognosis of CRC. In this study, based on the analysis of the TCGA data set, 484 differentially expressed IRGs were screened from 389 HCC and 39 normal tissues. By univariate regression analysis of differentially expressed IRGs, 30 genes were detected to be significantly correlated with OS. To further study the regulatory mechanisms of prognostic IRGs, a tumor-related TF-mediated network was established to reveal key TFs that can regulate these IRGs. Studies have shown that CBX7 played an important role in gastric and pancreatic cancer (, ). In recent years, studies have found that CBX7 was a component of polycomb repressive complex 1, maintaining the stem cell–like characteristics of gastric cancer cells by activating the AKT pathway and down-regulating p16 (). MYH11 (also known as SMMHC) encodes a smooth muscle myosin heavy chain, which plays a key role in smooth muscle contraction. The inversion of the MYH11 locus is one of the most common chromosomal aberrations in acute myeloid leukemia (). The MYH11 gene has a single-nucleotide repeat sequence (C8) in the coding sequence, which may be a mutation target for cancer that exhibits microsatellite instability (MSI). The study found that compared with the low microsatellite unstable group, the incidence of MYH11 frameshift mutation was higher in patients with high microsatellite-unstable (MSI) gastric cancer and CRC (). Among these major hub genes, the study of CXCL12 is more comprehensive. It has been reported that the CXCL12/CXCR4 axis is related to tumor progression, angiogenesis, metastasis, and survival (). Recent studies have found that the activation of LMO2 is essential for the development of T-cell acute lymphoblastic leukemia (T-ALL) leukemia (). The SP1PR1 gene plays a role in regulating tumors. Targeting the SphK1/S1P/S1PR1 axis with specific drugs can reduce tumor progression caused by key proinflammatory cytokines, macrophage infiltration, and obesity (). FOXP3 is one of the key TFs controlling the development and function of regulatory T cells. FOXP3 has been extensively studied in human tumors, which is closely related to tumor immunity, and its correlation with T cells in tumors has recently been reported (). The relationship between several other genes and CRC is still unclear. Among them, the role of MAF, A2M, EPAS1, and NR3C1 in CRC is worthy of further investigation. A study of breast cancer showed that the enhanced expression of MAF can mediate bone metastasis of breast cancer, which can be used as a risk index for bone metastasis in breast cancer patients (). The proteins encoded by A2M are protease inhibitors and cytokine transporters. It can inhibit a variety of proteases, as well as inflammatory cytokines, thereby destroying the inflammatory cascade. Xu’s team found that EPAS1 gene is dysregulated in non–small cell lung cancer, which encodes hypoxia-inducible factor 2α and plays an important role in the progression of non–small cell lung cancer (). It is known that EPAS1 is regulated by DNA methylation transcription in CRC (), but its role in CRC remains to be studied. NR3C1 encodes a glucocorticoid receptor, which can act both as a TF that binds to the glucocorticoid response element in the promoter of the glucocorticoid response gene to activate its transcription and as a regulator of other TFs. Further experimental evidence on the function of these genes in CRC may be of great help to our understanding of the progress of CRC.

In recent years, established a 20-gene prognosis model, which has a good predictive function for CRC prognosis. Another study also constructed a novel four-gene signature for CRC OS prediction based on gene expression data from TCGA, COAD, and READ data sets (). A recent study exploring the prognostic value of immune cells in the CRC tumor microenvironment determined that tumor-infiltrating immune cells is highly correlated with the progression and prognosis of CRC (). However, these studies do not fully explore the relationship between immune genes and the prognosis of CRC. Our study has the following advantages. First, we used a specialized immunological database to analyze as many IRGs as possible. To our knowledge, this is the first study to explore the relationship between a large number of IRGs and the prognosis of patients with CRC. Second, we obtained some immune-related prognostic genes and established a novel prognostic model related to immunity. This prognostic model showed excellent performance in the prediction of OS based on the TCGA database. According to the in-depth analysis, the immune-related prognostic model was demonstrated to be an independent prognostic indicator after adjusting for other clinical factors. These results indicated that the immune-related prognosis model can be used as an effective marker for the prognosis of CRC patients.

The characterization of immune infiltration is of great significance for studying the interaction between tumor and immunity. Therefore, we explored the relationship between immune-related prognostic models and immune cell infiltration to reflect the state of the immune microenvironment. According to the TIMER database, we found that high-risk patients had higher levels of CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells of infiltration. These results confirmed and extended the discovery that the heterogeneity of immune infiltration is important for the progression of CRC. A recent study reported that the colonic cancer microenvironment uses dendritic cells’ plasticity to support cancer progression by enhancing the release of the inflammatory chemokine CXCL1 (), which is consistent with our results. Neutrophils contribute to the activation, regulation, and effect of immune cells (). Existing research reported that tumor-associated neutrophils in CRC produce matrix metalloproteinase 9 vascular endothelial growth factor and hepatocyte growth factor to promote tumor invasion and angiogenesis. In addition, neutrophils also promote the spread of tumor cells by capturing tumor cells in the circulation, thereby promoting their migration to distant places (). Studies have reported that macrophages are associated with CRC progression (). Tumor-associated macrophages (TAMs) can induce EMT processes to enhance CRC migration, invasion, and circulating tumor cell (CTC)-mediated metastasis (). The immune model can indicate the infiltration of immune cells to some extent. It may be a promising way to cure CRC by broadening the relationship between immune cells and tumor progression.

Current research provides novel insights into the CRC immune microenvironment and immunotherapy. We conducted functional studies on selected genes to confirm their clinical value. However, the limitation of this study is that it is a retrospective study. Therefore, further prospective research is needed. On the one hand, the predictive capability of this model in CRC requires further testing with the goal of better prognostic stratification and treatment management. On the other hand, we need to further study the biological functions of the 18 IRGs through a series of experiments.

In short, through comprehensive analysis, many IRGs were found to be significantly related to the prognosis of CRC. Besides, we constructed a novel immune-related prognosis model as an independent prognostic indicator of CRC. This prognostic model can also indicate the infiltration of immune cells and prove its key role in the TIME. The current research has deepened our understanding of IRGs in CRC and provided new potential prognostic and therapeutic biomarkers.

Statements

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

YQ conceived the research. YQ and JW designed the research process, conducted bioinformatics analysis, downloaded and collated the data in the article, and wrote the first draft of the article. YQ, JW, and WL performed statistical analysis on the data. WL revised the article strictly to obtain the necessary knowledge and administrative support. FS, MH, KJ, DF, XZ, XK, QX, YH, and KD reviewed and edited the manuscript. All authors read and approved the final manuscript.

Funding

This study was funded by grants from the National Natural Science Foundation of China (81772545, 81802750, 81672916, and 81702331).

Acknowledgments

The authors would like to thank the TCGA, ImmPort, Cistrome Cancer, and TIMER databases for the availability of the data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.619611/full#supplementary-material

Supplementary Figure 1

The results of Kaplan–Meier analysis of 10 hub genes.

Supplementary Figure 2

Construction of an immune-related prognostic signature for CRC. (A) The risk score distribution of CRC patients in The Cancer Genome Atlas (TCGA) database. (B) Survival status and duration of patients. (C) Heatmap of the expression of 12 immune-related genes in CRC patients. (D) Survival curves for the low-risk and high-risk groups. (E) The receiver operating characteristic curve (ROC) analysis predicted overall survival using the risk score. The forecast time is 1, 3, and 5 years.

References

Summary

Keywords

colorectal cancer, immune-related gene, immune prognostic signature, TCGA, tumor immune microenvironment

Citation

Qian Y, Wei J, Lu W, Sun F, Hwang M, Jiang K, Fu D, Zhou X, Kong X, Zhu Y, Xiao Q, Hu Y and Ding K (2021) Prognostic Risk Model of Immune-Related Genes in Colorectal Cancer. Front. Genet. 12:619611. doi: 10.3389/fgene.2021.619611

Received

20 October 2020

Accepted

15 January 2021

Published

04 March 2021

Volume

12 - 2021

Edited by

Xiaoming Xing, The Affiliated Hospital of Qingdao University, China

Reviewed by

Edmund Ui-Hang Sim, Universiti Malaysia Sarawak, Malaysia; Xueqiu Lin, Stanford University, United States

Updates

Copyright

*Correspondence: Kefeng Ding,

These authors share first authorship

This article was submitted to Cancer Genetics, a section of the journal Frontiers in Genetics

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics