Identification of Molecular Subtypes and a Prognostic Signature Based on Inflammation-Related Genes in Colon Adenocarcinoma

Abstract

Both tumour-infiltrating immune cells and inflammation-related genes that can mediate immune infiltration contribute to the initiation and prognosis of patients with colon cancer. In this study, we developed a method to predict the survival outcomes among colon cancer patients and direct immunotherapy and chemotherapy. We obtained patient data from The Cancer Genome Atlas (TCGA) and captured inflammation-related genes from the GeneCards database. The package "ConsensusClusterPlus" was used to generate molecular subtypes based on inflammation-related genes obtained by differential expression analysis and univariate Cox analysis. A prognostic signature including four genes (PLCG2, TIMP1, BDNF and IL13) was also constructed and was an independent prognostic factor. Cluster 2 and higher risk scores meant worse overall survival and higher expression of human leukocyte antigen and immune checkpoints. Immune cell infiltration calculated by the estimate, CIBERSORT, TIMER, ssGSEA algorithms, tumour immune dysfunction and exclusion (TIDE), and tumour stemness indices (TSIs) were also compared on the basis of inflammation-related molecular subtypes and the risk signature. In addition, analyses of stratification, somatic mutation, nomogram construction, chemotherapeutic response prediction and small-molecule drug prediction were performed based on the risk signature. We finally used qRT-PCR to detect the expression levels of four genes in colon cancer cell lines and obtained results consistent with the prediction. Our findings demonstrated a four-gene prognostic signature that could be useful for prognostication in colon cancer patients and designing personalized treatments, which could provide new versions of personalized management for these patients.

Keywords: colon adenocarcinoma; immune infiltration; inflammation; molecule subtypes; signature.

Figure 1

(A) Volcano plot of 66 up-regulated and 54 down-regulated IRGs in COAD (FDR < 0.05 and |logFC| > 1). (B) Heatmap of 120 DE-IRGs between normal colon and COAD tissues. (C) The top ten enriched terms in GO analysis belonged to BP, CC, and MF for DE-IRGs. (D) The top thirty enriched terms in KEGG analysis. (E) The correlations between the top ten up-regulated and down-regulated IRGs. (F) PPI network of the DE-IRGs according to the STRING database. (G) The hub genes obtained from “cytohubba” plugin. (H, I) The two modules obtained from “MCODE” plugin. IRGs, Inflammation-related genes; COAD, Colon adenocarcinoma; FDR, False discovery rate; FC, Fold change; DE-IRGs, Differentially-expressed IRGs; GO, Gene Ontology; BP, Biological process; CC, Cell component; MF, Molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, Protein-protein interaction.

Figure 2

(A) Forest plot of ten prognostic-related DE-IRGs through univariate Cox analysis. (B) The correlations between the ten genes. (C) Consensus clustering matrix when k = 2. (D) Consensus clustering CDF with k valued 2 to 9. (E) Relative change in area under CDF curve for k = 2. (F) KM curve of the survival difference between cluster 1 and cluster 2. (G) Heatmap of the ten genes between the two clusters and the correlations of the clusters and clinical parameters. Immune cell infiltration using CIBERSORT (H), immune and stromal scores using ESTIMATE (I), the expression of MHC molecules (J), angiogenic activity, mesenchymal-EMT, tumourigenic cytokines and stemness scores (K), five common immunoinhibitors (L), and TIDE score (M) between the two clusters. CDF, Cumulative distribution function; KM, Kaplan–Meier; EMT, Epithelial-mesenchymal-transition; TIDE, Tumour Immune Dysfunction and Exclusion.

Figure 3

(A) Forest plot of the four genes selected in the signature through multivariate Cox analysis. (B) Coefficients of the four genes included in the signature. (C) The correlations between the signature and the four genes. (D) Heatmap of the association between the expression levels of the four genes and clinicopathological features. Survival analysis, heatmap, survival status accompanied with the risk score and ROC analysis in TCGA cohort (E) and GSE17538 cohort (F). The signature was an independent risk factor for COAD patients in TCGA cohort according to univariate (G) and multivariate Cox analysis (H). The differences of the risk score between different groups according to clinicopathological features, e.g., clusters (I), tumour stage (J), lymph node status (K), and metastasis (L). (M) Nomogram based on risk score and age. (N) Calibration plots of the nomogram for predicting the probability of 3.5-, 5- and 7.5-year survival. ROC, Receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Figure 4

(A) Pathways related to tumour development and progression enriched in the high-risk group. (B) JAK-STAT signaling pathway was the most relevant KEGG pathway in the high-risk group. (C) Multiply pathways associated with immune, chemokine and MHC molecules enriched in the high-risk group. Immune cell infiltration and immune-related functions or pathways (D), immune and stromal scores (E), immune cell infiltration using TIMER (F) and CIBERSORT (G), MHC molecules expression level (H), five common immunoinhibitors (I) and TIDE score (J) between the high- and low-risk groups. (*P < 0.05; **P < 0.01; ***P < 0.001; ns, not significant).

Figure 5

(A) Pathways related to angiogenesis, EMT, cytokine-cytokine receptor interaction and stemness enriched in the high-risk group. (B) Differences of angiogenic activity, mesenchymal-EMT, tumourigenic cytokines and stemness scores between the high- and low-risk groups. (C) The correlation of the risk score and angiogenic activity, mesenchymal-EMT, tumourigenic cytokines and stemness scores. (D) Differences of TSIs between the two groups. TSIs, Tumour stemness indices. (*P < 0.05; **P < 0.01; ***P < 0.001; ns, not significant).

Figure 6

Waterfall maps of the somatic mutations in the high-risk group (A) and the low-risk group (B). Heatmap of co-occurrence and mutually exclusive mutations of the differently mutated genes in the high-risk group (C) and the low-risk group (D). *p < 0.01. (E) Comparison of TMB between the high- and low-risk groups. (F) Difference in overall survival between high TMB and low TMB groups. (G) Difference in overall survival based on TMB and risk score. (H) Mutation rates of four genes (TIMP1, IL13, PLCG2, BDNF) in COAD patients from the cBioPortal database. (ns, not significant).

Figure 7

(A) The differences in the chemotherapy response of common chemotherapy drugs between the high- and low-risk groups. (B) Differentially expressed genes between the high- and low-risk groups. (C) The 3D structure of six potential target drugs screened out from the cMap database. (*P < 0.05; **P < 0.01; ***P < 0.001; NS, not significant).

Figure 8

The mRNA expression levels of TIMP1 (A), IL13 (B), BDNF (C) and PLCG2 (D) in different cell lines (NCM460, HCT116, SW480, HT29, LOVO, RKO, DLD-1) were measured by qRT-PCR. Results were normalized to reference gene GAPDH. Data are shown as the mean ± SEM, two-tailed unpaired t test was used for statistical calculation for each marker, n=3 independent experiments. (*P < 0.05, **P < 0.01, ***P < 0.001; ns, not significant).