Development of a prognostic index of colon adenocarcinoma based on immunogenomic landscape analysis.

BACKGROUND: Colon adenocarcinoma (COAD) is one of the most commonly diagnosed cancers, and it is closely related to the immune microenvironment. Considering that immunotherapy is not effective for all COAD patients, it is necessary to identify the effective population before administering treatment. In this study, we established an independent prognostic index based on immune-related genes (IRGs), in order to evaluate the clinical outcome of COAD.
METHODS: The gene expression profiles and IRGs taken from The Cancer Genome Atlas (TCGA) and Immunology Database and Analysis Portal (ImmPort), respectively, were integrated in order to identify the differentially expressed IRGs. Functional enrichment analysis was conducted and the prognostic value of survival-related IRGs was determined. Based on Cox regression analysis, the IRG-based prognostic index (IRGPI) was established, and the model was evaluated and applied.
RESULTS: A total of 51 differentially expressed survival-related IRGs were identified. The most significant signaling pathway was "cytokine-cytokine receptor interaction". The index established herein was based on 12 survival-related IRGs, and it was highly accurate in monitoring prognosis. Moreover, the IRGPI was significantly correlated with multiple clinicopathologic factors, as well as with the infiltration of immune cells.
CONCLUSIONS: An independent IRGPI was established in order to assess the immune status and tumor prognosis in COAD patients. This index can serve as a robust biomarker in clinical prognosis applications, including cancer immunotherapy. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

According to the Global Cancer Report released at the end of 2018, colon adenocarcinoma (COAD) is one of the most commonly diagnosed cancers, and one of the top three malignant tumors in terms of morbidity and mortality (1). Considering that the early symptoms are not visible, most patients diagnosed with COAD are in the middle or late stages of the disease, which impedes treatment and timely surgical intervention. Therefore, despite the increasing number of targeted drugs (2,3), the 5-year average survival rate of COAD patients in China is only 57%. To improve this rate, the prognosis status of patients should be monitored closely, and the treatment plans should be timely adjusted. This necessitates the development of an independent prognostic index that allows for the evaluation of COAD clinical outcomes.

Evidence of the close relation between the immune microenvironment and tumor development has substantially increased over the past few years (4-7). Consequently, tumor immunotherapy has emerged as a promising technique of treatment, and it has garnered the attention of many researchers (8,9). This technique relies on the activation of immune cells that can trigger anti-tumor immune responses, leading to the elimination of tumor residual lesions, the inhibition of tumor growth, and the breaking of immune tolerance (10). Unfortunately, immunotherapy is not effective for all tumor patients. Therefore, clinical research efforts are currently focused on the identification of patients who can be efficiently treated using immunotherapy techniques. This may be accomplished by determining the type and proportion of immune cells and factors in the tumor immune microenvironment, which are in turn controlled by gene expression. The recent progress in gene sequencing technology has allowed for the design of new and sensitive gene-based biomarkers (11,12) that may be used to improve the prognosis of COAD. The comprehensive analysis of the relevance between the immune-related genes (IRGs) and survival is conducive to elucidating the potential prognostic value of these genes. For example, an immuno-genomic study of tumor immune microenvironment has previously shown that mRNA-seq is effective in predicting clinical response (5).

In this study, we identified the differentially expressed IRGs in COAD patients by relating the gene expression profiles to the IRGs from different online databases. The functional enrichment analysis and prognostic value of the survival-related IRGs were also determined. The prognostic index was finally established based on Cox regression analysis, followed by model evaluation and application. In theory, the robust prognostic biomarker reported herein can be used to predict the clinical outcomes of therapeutic methods, including cancer immunotherapy.

Methods

Acquisition of gene expression data and IRGs

The RNA-seq data of COAD samples was taken from The Cancer Genome Atlas (TCGA) database (13) that is published in the Genomic Data Commons (GDC) portal (https://portal.gdc.cancer.gov/) (14), and the clinical information of these samples was downloaded and extracted by Perl. The collected dataset is comprised of 473 COAD samples and 41 adjacent normal samples. In addition, a comprehensive list of IRGs (1,811 in total) was downloaded from the Immunology Database and Analysis Portal (ImmPort) (https://www.immport.org/home) (15). These genes were divided into 17 categories based on their molecular functions (MFs) (e.g., cytokine, interleukin, and natural killer cell cytotoxicity).

Identification of differentially expressed genes (DEGs)

The DEGs were identified by comparing the COAD and adjacent normal samples using R package “limma” (16). Specifically, the genes exhibiting |log2(fold change)| >1 and false discovery rate (FDR) <0.05 were considered as DEGs. Afterwards, the differentially expressed IRGs were extracted from the list of identified DEGs. Heat maps and volcano plots of DEGs and differentially expressed IRGs were generated by R package pheatmap and ggplot2, respectively (17,18).

Screening of survival-related IRGs and evaluation of their prognostic value

A survival analysis of all of the differentially expressed IRGs was performed using the R package “survival” (19), and the survival-related IRGs were determined by univariate Cox analysis (P<0.05). In order to explore the potential molecular mechanisms of differentially expressed IRGs and survival-related IRGs, functional enrichment analysis were conducted on the basis of Kyoto encyclopedia of genes and genomes (KEGG) and gene ontology (GO) database, including biological process (BP), MF and cellular component (CC) (20,21), using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) 6.8 (22). The GO terms and KEGG signaling pathways were enriched with an adjusted P value of <0.05. Moreover, the top 10 KEGG signaling pathways were drawn into a bubble chart using the R package “ggplot2” (18). To determine the prognostic value of survival-related IRGs, a forest plot of these genes was constructed using the hazard ratio (HR) as an indicator, based on univariate Cox analysis.

Development of the IRG-based prognostic index (IRGPI)

Survival-related IRGs were submitted for multivariate analysis to develop the IRGPI, and integrated IRGs was kept as independent prognostic indicators. Specifically, the IRGPI was established by multiplying the expression values with the Cox regression coefficient (23). It should be noted that the survival-related IRGs included in the IRGPI are referred to as hub IRGs.

Assessment of IRGPI and genetic alteration analysis

Patients were divided into high- and low-risk groups, depending on whether their IRGPI values are greater or less than the median value, respectively. The corresponding Kaplan-Meier (K-M) survival curves were constructed to demonstrate the overall survival (OS) within different risk groups. The ROC curve was also established in order to assess the sensitivity and specificity of the model. Additionally, univariate and multivariate analyses of survival were carried out for both IRGPI and clinicopathologic factors. These analyses were performed using the R package “survival” (19). We also explored the correlation between the expression of hub IRGs and the clinicopathologic factors. Moreover, genetic alteration analysis of hub IRGs was conducted on cBioPortal (https://www.cbioportal.org/) (24).

Evaluation of immune cell infiltration in tumors

The Tumor Immune Estimation Resource (TIMER) was used to systematically analyze immune cell infiltration in 23 types of cancer listed in the TCGA (25). Statistical methods were adopted to estimate the abundances of six immune cell infiltrates, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. The obtained results were validated by pathological estimations. Subsequently, the immune cell infiltration levels in COAD samples were downloaded, and the correlation between IRGPI and immune cell infiltration was calculated based on Pearson correlation analysis.

Construction of the transcription factor (TF)-mediated regulatory network

TFs play an essential role in controlling gene expressions, including IRG, and thus, they can affect the immune function of the body. Consequently, it is necessary to identify the particular TFs that can potentially regulate hub IRGs. For this purpose, the list of TFs was first downloaded from the Cistrome Cancer database (http://cistrome.org/CistromeCancer/) (26). Then, the differentially expressed TFs were extracted from all DEGs. Subsequently, we carried out the correlation analysis of differentially expressed TFs and hub IRGs using the R software. The correlativity was considered reliable if |correlation value| >0.6 and P<0.05. As a result, we submitted them to build a TF-mediated regulatory network via Cytoscape software. To screen for the most important TF in this network, topological analysis was performed using the Network Analyzer plug-in (27).

Statistical analysis

To evaluate IRGPI effectiveness, the R package “survival” was used to calculate the area under curve (AUC) of the receiver operating characteristic (ROC) curve (28), and the R package “ggplot2” was used to generate boxplots (18). An independent t-test was conducted to test the differences among diverse clinical parameters. These differences were considered to be statistically significant if P<0.05.

Results

Acquisition of differentially expressed IRGs

In order to identify the DEGs, the genomic data of 473 COAD and 41 adjacent normal samples were compared using R package “limma”. Based on the cut-off criteria of |log2(fold change)| >1 and FDR <0.05, a total of 6,477 DEGs were detected, including 4,561 up-regulated genes and 1,916 down-regulated genes (). With the list of IRGs, 467 differentially expressed IRGs were extracted from all DEGs, including 179 up-regulated IRGs and 288 down-regulated IRGs ().

Figure 1

Differentially expressed IRGs. The heatmap (A) and volcano plot (B) of genes that are differentially expressed in COAD samples, compared to adjacent normal samples; red and blue dots represent up-regulated and down-regulated DEGs, respectively; the heatmap (C) and volcano plot (D) of differentially expressed IRGs. Red and blue dots represent up-regulated and down-regulated differentially expressed IRGs, respectively. IRGs, immune-related genes; COAD, colon adenocarcinoma; DEGs, differentially expressed genes; FDR, false discovery rate.

Functional enrichment analysis of differentially expressed IRGs

To explore the biological functions of differentially expressed IRGs, functional enrichment analysis was conducted using the DAVID 6.8. The results illustrated in show that “immune response”, “extracellular region”, and “antigen binding” are the most relevant BP, CC, and MF of genes, respectively. Furthermore, the KEGG enrichment analysis presented in indicates that the differentially expressed IRGs are mostly correlated with “cytokine-cytokine receptor interaction”.

Figure 2

Functional enrichment analysis of differentially expressed IRGs. (A) GO enrichment analysis. Blue, red, and green parts represent BP, CC, and MF, respectively; (B) the top 10 KEGG signaling pathways. IRGs, immune-related genes; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto encyclopedia of genes and genomes.

Functional enrichment analysis of survival-related IRGs

As key prognostic indicators in applications of clinical tumor surveillance, molecular biomarkers related to IRGs could be used to evaluate the therapeutic efficiency of cancer immunotherapy. Based on univariate Cox analysis, 51 survival-related IRGs were distinguished (P<0.05). As shown in , the most relevant BP, CC, and MF terms of these genes are “regulation of response to stimulus”, “extracellular space”, and “receptor ligand activity”, respectively. Again, “cytokine-cytokine receptor interaction” was found to be the most significant signaling pathway ().

Figure 3

Functional enrichment analysis of survival-related IRGs. (A) GO enrichment analysis. Blue, red, and green parts represent BP, CC, and MF, respectively; (B) the top 10 KEGG signaling pathways. IRGs, immune-related genes; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto encyclopedia of genes and genomes.

Prognostic values of survival-related IRGs

In order to better define the characteristics of survival-related IRGs, the prognostic values of these genes were calculated. The resulting forest plot presented in demonstrates that 40 genes are characterized by HR >1, whereas 11 genes have HR <1. This suggests that most survival-related IRGs were risk factors affecting the prognosis of COAD.

Figure 4

Prognostic values of survival-related IRGs. P values <0.05 are considered to be statistically significant. IRGs, immune-related genes.

Development of IRGPI

In view of the results of the multivariate Cox regression analysis, the IRGPI was developed. This index was based on the following formula: risk score = [expression level of CD1B × (–2.09883)] + [expression level of SLC10A2 × 0.718209] + [expression level of FGF2 × 0.354982] + [expression level of CCL28 × (–0.1162)] + [expression level of IGHV4-31 × 0.00944] + [expression level of IGKV1-6 × 0.006782] + [expression level of ESM1 × 0.199876] + [expression level of TNFSF12 × 0.09248] + [expression level of UCN × 0.376204] + [expression level of VIP × 0.066404] + [expression level of GLP2R × (–4.08525)] + [expression level of IL1RL2 × 0.156916]. Depending on the scores and on the survival time, the patients were divided into high- (IRGPI > median value) and low-risk (IRGPI < median value) groups. As shown in , the high-risk group exhibits more deaths than the low-risk one.

Figure 5

Development of the IRGPI. (A) Risk scores of patients in different groups; (B) survival status of patients in different groups; (C) heat map of hub IRG expression profiles. IRGPI, immune-related genes-based prognostic index.

Evaluation of IRGPI

The calculated IRGPI values were used to determine the survival probability of COAD patients. As shown in , the OS of high-risk group patients was found to be less than that of the low-risk group patients (P<0.001). To assess the predictive accuracy of the IRGPI established herein, the ROC curve was generated, and the AUC was calculated (). This area was found to be 0.855, which indicates that the proposed IRGPI has excellent potential in monitoring prognosis. Moreover, the results of both univariate and multivariate analyses () demonstrate that IRGPI is significantly correlated with survival, and thus, it can be used as an independent indicator when multiple clinicopathologic factors are taken into account.

Figure 6

The evaluation of the IRGPI. (A) The K-M survival analysis of the IRGPI. Patients with high-risk scores showed poor OS, compared to those with low-risk scores; (B) survival-dependent ROC curve validation of the prognostic value of the IRGPI. The X-axis represents the false positive rate, while Y-axis represents the true positive rate. IRGPI, immune-related genes-based prognostic index; K-M, Kaplan-Meier; OS, overall survival; ROC, receiver operating characteristic; AUC, area under curve.

Table 1

Univariate and multivariate analysis of COAD

Variables	Univariate analysis			Multivariate analysis
	HR (95% CI)	P value		HR (95% CI)	P value
Age	1.024 (0.994–1.055)	0.117		1.029 (0.996–1.063)	0.090
Gender (male/female)	1.020 (0.524–1.987)	0.953		1.102 (0.528–2.299)	0.795
Tumor status (with tumor/tumor free)	2.673 (1.153–6.198)	0.022		1.662 (0.673–4.105)	0.271
Pathological stage	2.686 (1.804–3.998)	<0.001		0.727 (0.197–2.687)	0.633
T stage	3.637 (1.890–7.000)	<0.001		2.612 (1.157–5.895)	0.021
M stage	6.103 (3.110–11.976)	<0.001		5.443 (0.858–34.516)	0.072
N stage	2.143 (1.430–3.211)	<0.001		1.408 (0.689–2.879)	0.348
IRGPI	1.339 (1.226–1.463)	<0.001		1.354 (1.212–1.513)	<0.001

COAD, colon adenocarcinoma; HR, hazard ratio; IRGPI, immune-related genes-based prognostic index.

To further assess the clinical value of IRGPI, the relationships between the hub IRGs implicated in IRGPI and the clinicopathologic factors were analyzed. As shown in the results of , tumor status and M stage involved the expression of more genes. Besides, as an independent indicator, IRGPI had statistically significant difference in pathological stage (), T stage (), N stage (), and M stage (), which suggested that IRGPI could make an accurate prediction in different pathological stages of COAD tumors. However, IRGPI showed no significant difference in age (), gender () and tumor status ().

Table 2

The correlation between the expression of the IRGs and the clinical factors in COAD

Genes	Age (>65/≤65)			Gender (male/female)			Tumor status (with tumor/tumor free)			Pathological stage (III&IV/I&II)			T stage (T3–4/T1–2)			M stage (M1/M0)			N stage (N1–3/N0)
	t	P		t	P		t	P		t	P		t	P		t	P		t	P
CD1B	–0.155	0.877		1.717	0.087		1.519	0.13		2.849	0.005		1.579	0.118		4.658	8.85E–06		2.653	0.008
SLC10A2	0.775	0.439		1.711	0.089		0.859	0.391		–0.193	0.847		–0.871	0.385		2.211	0.028		–0.256	0.798
FGF2	–1.014	0.311		–1.2	0.231		4.287	2.98E–05		0.449	0.654		0.354	0.725		2.743	0.007		0.152	0.879
CCL28	–2.463	0.014		0.09	0.928		–2.198	0.029		2.192	0.029		0.612	0.542		3.013	0.003		2.167	0.031
IGHV4-31	–1.508	0.133		–2.096	0.037		–2.717	0.007		1.675	0.095		0.573	0.568		2.366	0.02		1.451	0.148
IGKV1-6	–0.071	0.944		–0.478	0.633		–0.706	0.481		0.459	0.646		–0.899	0.37		1.955	0.054		0.353	0.725
ESM1	1.269	0.205		0.033	0.973		2.314	0.021		–1.638	0.103		–1.866	0.065		–1.132	0.263		–1.843	0.067
TNFSF12	1.917	0.056		0.864	0.388		1.455	0.147		–0.761	0.447		–1.383	0.169		0.469	0.641		–1.002	0.317
UCN	–1.765	0.079		0.298	0.766		–5.221	3.80E–07		–2.052	0.041		–0.941	0.349		–1.135	0.262		–1.866	0.064
VIP	0.993	0.321		0.907	0.365		2.48	0.014		–1.672	0.096		–2.215	0.029		–0.586	0.56		–2.071	0.04
GLP2R	–0.88	0.379		0.408	0.684		–0.402	0.688		–1.293	0.197		0.363	0.718		0.318	0.751		–1.323	0.187
IL1RL2	–0.956	0.34		–0.078	0.938		0.722	0.471		–1.045	0.297		–0.257	0.798		–1.622	0.111		–0.481	0.631
IRGPI	1.195	0.233		–0.289	0.773		1.129	0.26		–2.433	0.016		–5.25	3.53E–07		–2.073	0.043		–2.457	0.015

IRGs, immune-related genes; COAD, colon adenocarcinoma; IRGPI, immune-related genes-based prognostic index.

Figure 7

The relationships between the prognostic index and (A) pathological stage, (B) T stage, (C) N stage, (D) M stage, (E) age, (F) gender, and (G) tumor status.

Mutation landscape of hub IRGs included in IRGPI

Considering the significant clinical value of hub IRGs, their molecular characteristics were comprehensively explored. Specifically, genetic alteration analysis of these genes was carried out using the cBioPortal. The obtained results indicate that the most frequent type of genetic alteration is the elevated expression of mRNA (). For most IRGs (75%), the genetic alteration rate was found to be ≥5%, with the ESM1 and IL1RL2 genes showing the highest alteration frequencies ().

Figure 8

Genetic alterations associated with hub IRGs. (A) The total alteration frequency of hub IRGs; (B) the landscape of genetic alterations of hub IRGs in COAD patient. IRGs, immune-related genes; COAD, colon adenocarcinoma.

Clinical application of IRGPI

In order to estimate whether the IRGPI could accurately reflect the status of tumor immune microenvironments, we analyzed the relationship between IRGPI and the abundance of immune cell infiltration. The obtained results demonstrate that IRGPI is significantly correlated with 5 types of immune cells, including B cells (), CD4 T cells (), dendritic cells (), macrophages (), and neutrophils (). Moreover, there was no significant correlation between IRGPI and CD8 T cells ().

Figure 9

Relationships between the prognostic index and the abundance of immune cell infiltration. (A) B cells, (B) CD4 T cells, (C) dendritic cells, (D) macrophages, (E) neutrophils, and (F) CD8 T cells.

TF-mediated regulatory network

In order to explore the molecular mechanisms of hub IRGs, we studied their regulatory mechanism. Knowing that TFs play a key role in regulating the molecular network, the expression profiles of 318 TFs were examined, and 68 differentially expressed TFs were identified (). As a result, a TF-mediated regulatory network was constructed based on 15 TFs and 8 hub IRGs (). Afterwards, we performed a topological analysis on these TFs based on three centrality algorithms, including degree centrality, betweenness centrality and closeness centrality. As shown in , MAF, MYH11, and NR3C1 ranked the top 5 TFs in all three centrality analyses, and thus, they were considered as key TFs in the TF-mediated regulatory network.

Figure 10

TF-mediated regulatory network. The heatmap (A) and volcano plot (B) of differentially expressed TFs. Red and blue dots represent up-regulated and down-regulated differentially expressed TFs, respectively; (C) the regulatory network. Green triangles represent TFs, whereas red and blue circles represent high-risk and low-risk IRGs, respectively. IRGs, immune-related genes; TF, transcription factor.

Table 3

The topological analysis of TFs in regulating the hub IRGs

TF name	Degree centrality	TF name	Betweenness centrality	TF name	Closeness centrality
MAF	4	MAF	0.202086	KLF4	1
MYH11	4	MYH11	0.129743	MAF	0.487805
NR3C1	4	FOXP3	0.114229	MYH11	0.465116
CBX7	3	IRF4	0.1	FOXP3	0.444444
FOXP3	3	NR3C1	0.094054	NR3C1	0.444444
EPAS1	2	CBX7	0.079211	CBX7	0.444444
IRF4	2	EPAS1	0.038571	EPAS1	0.408163
BHLHE40	1	BHLHE40	0	BRCA1	0.377358
BRCA1	1	BRCA1	0	CENPA	0.377358
CENPA	1	CENPA	0	EZH2	0.377358
EZH2	1	EZH2	0	LMO2	0.377358
KLF4	1	KLF4	0	NCAPG	0.377358
LMO2	1	LMO2	0	BHLHE40	0.307692
NCAPG	1	NCAPG	0	IRF4	0.28169
PBX1	1	PBX1	0	PBX1	0.266667

TF, transcription factor; IRGs, immune-related genes.

Discussion

The development of targeted drugs and cancer immunotherapy has shifted COAD treatment from surgery and chemoradiotherapy to precise and individualized therapies (29). However, due to the limitations of immunotherapy treatments, it is necessary to identify the effective population and predict the clinical outcome beforehand. This may be accomplished by analyzing the immune genes and prognostic status of patients. In this study, we established an IRGPI that may be used to evaluate immune response in COAD patients, and thus, serves as a new biomarker for cancer immunotherapy.

Our results demonstrate that the most relevant BP of differentially expressed IRGs is “immune response”, and that the survival-related IRGs are highly associated with “the regulation of response to stimulus”. Based on KEGG enrichment analysis, the survival-related IRGs were significantly correlated with “cytokine-cytokine receptor interaction”, which suggested that the response of the immune system to stimuli is related to the expression of IRGs and patients’ survival. Furthermore, assessments of the prognostic efficiency of survival-related IRGs showed that most genes are characterized by high HR values, and thus, they affect the prognosis of COAD patients. This explains the variations in the efficacy of the same treatment between different patients.

Currently, the methods used to identify potential mechanisms and prognostic biomarkers are based on different types of bioinformatic analyses. For example, Dalerba et al. used the bioinformatics approach in combination with clinical-grade diagnostic assays to show that a subgroup of patients benefits from adjuvant chemotherapy (30). Other researchers used gene expression signatures based on mRNA, miRNA, or lncRNA expression profiles to predict the prognosis of COAD (31-33). In this study, we propose a biomarker based on IRGs, which is the most suitable way to reflect the immune status and tumor prognosis in patients. To establish this biomarker, the expression levels of survival-related IRGs were integrated using Cox regression analysis. Ultimately, 12 of 51 extracted genes were used to construct the predictive model that was used to calculate the risk scores of COAD patients. Based on these scores, the patients in the analyzed sample were divided into high- and low-risk groups.

The model was also used to establish an IRGPI. The reliability and efficiency of this index in terms of survival prediction was confirmed by ROC curve analysis. Univariate and multivariate analyses further prove the significant correlation between IRGPI and survival. Although the univariate analysis yields 5 independent indicators, only one (IRGPI) could be used in multivariate analysis. To test the feasibility of the IRGPI, the relationship between hub IRGs and multiple clinicopathological factors was also analyzed. The results indicate that as a comprehensive index, IRGPI exhibits stronger correlation with multiple pathological stages than single genes. Furthermore, the index is capable of accurate survival prediction in the pathological, T, N, and M stages.

Considering that cancer and immunity are closely associated, clinical treatments must take tumor immune microenvironments into account. Assessments of the correlation between the IRGPI and immune cell infiltration can be used to reflect the status of such microenvironments. Herein, we show that the IRGPI is significantly correlated with the infiltration levels of B cells, CD4+ T cells, dendritic cells, macrophages, and neutrophils (P<0.05). These immune cells play a very important role in COAD prognosis. For example, Gu et al. have recently shown that the elimination of pathogenic CD4+ T-cells and the induction of antitumor CD8+ T-cell activity can suppress colon carcinogenesis (34). Also, Vo et al. demonstrate that dendritic cell and lenalidomide combination therapy can effectively enhance antitumor immunity in a mouse COAD model (35). Moreover, according to Vinnakota et al., M2-like macrophages induce COAD cell invasion via matrix metalloproteinases (36).

To explore the molecular mechanisms underlying hub IRGs activity, a TF-mediated regulatory network was established. This network is comprised of several TFs that can simultaneously regulate the expressions of multiple IRGs, such as MAF, MYH11, and NR3C1. The available research indicates that c-MAF is capable of controlling immune responses by regulating disease-specific gene networks and repressing IL-2 in CD4+ T cells (37). As for MYH11, it exhibits significantly different expressions in tumor samples, compared to normal tissues, and thus, it is considered as a potential tumor biomarker (38). Finally, the NR3C1 gene expressed during CD8+ T cell differentiation regulates the formation of memory-precursor cell fates (39). Other studies have also investigated the role of TFs in COAD prognosis. For instance, Wei et al. established a regulatory network of TF-miRNA-target genes to analyze disease progression and optimize treatment strategies (40). However, Mullany et al. showed that TF expression and related miRNAs influence the survival of patients diagnosed with colorectal cancer (41). Based on these studies and on our own results, we believe that TFs may play an important role in regulating the IRG expressions in COAD patients. Further investigation is needed to confirm this hypothesis.

Considering that the IRGPI established herein may be applied to accurately assess the immune status and monitor COAD prognosis in patients, it is actually an outstanding biomarker of the disease. This robust biomarker can be used in clinical prognosis applications, including cancer immunotherapy.