This article is mentioned in:

Introduction

Colorectal cancer (CRC) is the third most frequently diagnosed tumor, with the second highest mortality rate globally in 2018 (1). Due to its distinct molecular characteristics and genetic heterogeneity, patients with CRC often have different natural processes and clinical outcomes. Therefore, there is an urgent requirement to explore the molecular subtypes of CRC for precise therapy. Recently, several studies have reported some distinct subtypes of CRC. However, the subclasses generated from these studies either lack in the analysis of different molecular mechanisms (2,3) or have insufficient independent samples for validation (4,5). Thus, a more effective and credible subgroup is required for the treatment of patients with CRC.

In recent years, owing to the further understanding of the interactions between the immune system and tumor cell environment, immunotherapies such as anti-programmed cell death protein 1 (PD-1), anti-programmed cell death 1 ligand 1 (PD-L1) and anti-cytotoxic T-lymphocyte protein 4 (CTLA-4) have markedly revolutionized the therapeutic approaches in several cancer types, including lung cancer (6,7), melanoma (8,9), head and neck cancer (10), bladder cancer (11), kidney cancer (12) and CRC with deficient mismatch repair (13,14). An ~50% elevated response rate was found, which is a noteworthy improvement in patients who received ICI treatment due to the reversion of an immunosuppressive microenvironment (15,16). However, only a small number of CRC patients exhibit a clinical response to ICI treatment, with the majority of patients not benefiting from the therapy (17). The therapeutic effect of ICI agents is ascribed to several important factors, such as expression of PD-L1, proportion of tumor-infiltrating lymphocytes (TILs) and tumor mutation burden (TMB). High PD-L1 expression is an essential factor in the use of pembrolizumab, which is approved by the Food and Drug Administration (FDA) for patients with non-small cell lung cancer (NSCLC) (18). Several studies have also reported the vital roles of TIL and TMB in the response to ICI agents (19–21). In summary, the exploration of the CRC subtype that is suitable to receive ICI therapy is important and urgently required.

In the present study, immune-related gene expression profiles of four independent CRC cohorts were used to identify the potential subtype which may be sensitive to immune checkpoint blockade. Findings of the present study have implications for guiding clinical immunotherapy for patients with CRC.

Materials and methods

Acquisition of genomic data

Gene expression data of 662 CRC samples with overall survival time and status from TCGA cohort (https://gdc.cancer.gov) and 619 samples of three datasets (GSE103479, GSE38832 and GSE87211) (22–24) from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) were retrospectively acquired. TCGA cohort (mean age, 66.33; age range: 31–90; sex distribution: 355 male and 307 female patients); GSE103479 dataset (mean age: 64.26; age range: 35–88; sex distribution: 83 male and 72 female patients); GSE8832 dataset (mean age: 59.63; age range: 32–89; sex distribution: 69 male and 53 female patients) and GSE87211 dataset (mean age: 63.16; age range: 38–88; sex distribution: 182 male and 160 female patients). Clinical characteristics of these four cohorts are presented in Table SI. Finally, a total of 1,281 CRC samples with gene expression data and survival information from four independent datasets were obtained (Table SII). Complete clinical data were available for all patients included in the present study. Gene expression profiles were consistently normalized to reduce variance. Expression levels of genes with multiple probes were presented as the mean expression of all probes.

Immune-related prognostic signature used for molecular subtyping

A prognostic signature was established based on a previously conducted immune-related genes study (25). The aforementioned previous study comprehensively described the immune landscape of >10,000 samples, comprising 33 different cancer types, and integrated 160 immune-related signatures containing 2,995 immune genes. Univariate Cox proportional hazards model was used to evaluate the association between the expression of these genes and the overall survival of patients with CRC. Genes with P<0.01 were included in the signature to perform molecular subtyping.

NMF clustering analysis

Molecular subtyping was performed applying non-negative matrix factorization (NMF) (26,27). A binary matrix (A) describing the expression of immune-related genes (rows) across CRC samples (columns) was established. Subsequently, the expression matrix A was factorized into two non-negative matrices (W and H; A≈WH). Matrix H was applied to cluster samples into distinct subtypes. The values of cophenetic, dispersion and silhouette coefficients were used to select the optimal number of subtypes. The NMF clustering analysis was performed with the R package ‘NMF’ v.0.21.0 (28).

Prediction of response to ICI treatment

The tumor immune dysfunction and exclusion (TIDE) algorithm (29) was applied to predict potential distinct responses to ICI therapy. TIDE is a gene expression biomarker for predicting the response to immune checkpoint blockade in patients. A low TIDE prediction score represents weak potential immune escape, and therefore these patients would potentially exhibit a greater immune therapy response. TIL proportion was analyzed using the CIBERSORT algorithm (30), which is a useful analytical tool to provide an evaluation of the abundances of 22 immune cell types in a mixed cell population, using gene expression data.

Gene set enrichment analysis (GSEA)

Patients with CRC were partitioned into two groups, according to the subtyping results clustered with the molecular expression features. ‘DEseq2’ v.1.26.0 (31) and ‘limma’ v.3.34.8 (32) packages were applied to calculate the differential t statistics of the RNA sequencing and microarray data. The t statistic was used as the input to R function in the ‘fgsea’ v.1.12.0 package (http://bioconductor.org/packages/release/bioc/html/fgsea.html) to perform GSEA. The pathway annotation signatures from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases in Molecular Signatures Database (MSigDB) (33) were used.

Normalization of gene expression data

In the present study, the molecular subtyping of CRC samples was conducted using the NMF algorithm. The essential condition for performing the NMF approach was non-negative values, therefore, the gene expression data was compressed into the range from 0–1 in all distinct platforms to achieve data normalization.

Statistical analyses

Statistical analyses were conducted with R software 3.6.1 (https://cran.r-project.org). The difference in clinical characteristics between two subtypes were compared using the χ2 test, and differences of TIDE score, PD-L1 expression and TMB were compared using the Wilcoxon rank-sum test. Survivals plot were drawn using the Kaplan-Meier method and log-rank test for comparison. The association between CRC subtypes and prognosis was analyzed with univariate and multivariate Cox proportional hazards model in the R ‘survival’ package (v.2.41–3) (https://gitub.com/therneau/survival). P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of two CRC subtypes with distinct prognoses

Using the univariate Cox proportional hazard model, the association between the expression of 2,995 immune-related genes and the prognosis of patients with CRC in TCGA cohort was determined. Finally, 53 genes with P<0.01 (Table SIII) were included in the prognostic signature for molecular subtyping. NMF unsupervised clustering analysis of these 53 immune genes was performed and the results showed that the cophenetic, dispersion and silhouette coefficients harbored the maximum value for the factorization rank of two subtypes (Figs. 1A and S1). The survival analysis of these two subtypes indicated that the overall survival rate of patients in the high-risk subtype was lower compared with that of patients in the low-risk subtype (HR, 1.93; 95% CI, 1.32–2.87; log-rank P<0.001; Fig. 1B). Multivariate Cox regression analysis adjusted for age, gender, stage and histology remained statistically significant (HR, 1.56; 95% CI, 1.08–2.27; P=0.018; Fig. 1C).

Figure 1. Identification of two distinct subtypes of CRC with significant prognostic differences. (A) Association of cophenetic, dispersion and silhouette coefficients with the number of subtypes. (B) Kaplan-Meier plot of the two identified CRC subtypes (low-risk vs. high-risk). (C) Forest plot of multivariate Cox regression model of low-risk vs. high-risk subtypes with confounding clinical factors taken into consideration. CRC, colorectal cancer; COAD, colon cancer; READ, rectal cancer.

In order to verify whether the aforementioned subtyping approach would discover the same subtypes in other CRC cohorts, the same clustering approach was applied to analyze 3 independent CRC cohorts, including GSE103479 (n=155), GSE38832 (n=122) and GSE87211 (n=342). Consistently, two different subtypes were identified in all 3 additional cohorts and the high-risk subtype also exhibited a significantly poorer prognosis [GSE103479: HR, 2.09 (95% CI, 1.47–3.53; log-rank P=0.005); GSE38832: HR, 1.78 (95% CI, 1.25–3.16; log-rank P<0.001) and GSE87211: HR, 1.78 (95% CI, 1.25–3.16; log-rank P<0.001); Fig. 2].

Figure 2. Kaplan-Meier survival plots of low-risk and high-risk colorectal cancer subtypes in (A) GSE103479, (B) GSE38832 and (C) GSE87211 datasets.

Prediction of response to ICB therapy in the two CRC subtypes

Since the two CRC subtypes were identified by the application of 53 immune-related genes, it was hypothesized that these two subtypes may have distinct responses to ICB therapy. Thus, the TIDE algorithm was used to predict the potentially different responses to ICB therapy. In TCGA cohort, the TIDE prediction score was significantly lower in the low-risk subtype compared with that in the high-risk group (P<0.001; Fig. 3A), which was further validated in the three other independent cohorts (GSE103479, GSE38832 and GSE87211; all P<0.01; Fig. 3B-D). These findings revealed that the patients of the low-risk subtype may be more sensitive to ICB treatment.

Figure 3. Differences in TIDE prediction score between the two CRC subtypes in (A) The Cancer Genome Atlas and (B-D) three validation CRC cohorts. CRC, colorectal cancer; TIDE, tumor immune dysfunction and exclusion.

Differences in ICI response markers between the two CRC subtypes

PD-L1 expression, TIL and TMB are highly important clinical factors, due to their vital roles in immunotherapy responses. Therefore, the differences in these 3 factors were compared between the two CRC subtypes. PD-L1 expression was significantly higher in the low-risk subtype compared with that in the high-risk subtype in TCGA cohort (P<0.001), which was further validated in three other independent cohorts (GSE103479, GSE38832 and GSE87211; all P<0.05; Fig. 4A).

Figure 4. Differences in PD-L1 expression, TIL proportion and TMB in low-risk vs. high-risk CRC subtypes. Differences in (A) PD-L1 expression and (B) TIL proportion in TCGA and three independent validation cohorts, stratified by the two CRC subtypes. (C and D) Association of TMB with low-risk and high-risk CRC subtypes. CRC, colorectal cancer; TCGA, The Cancer Genome Atlas; PD-L1, programmed death-ligand 1; TIL, tumor-infiltrating lymphocyte; TMB, tumor mutation burden.

Results showed that the TIL proportion was significantly higher in the patients of the low-risk subtype compared with that in the high-risk group in TCGA cohort (P=0.003). This association was also corroborated in the three additional validation cohorts (GSE103479, GSE38832 and GSE87211; all P<0.05; Fig. 4B). The differences in the infiltration of 22 immune cells between the two identified CRC subtypes in the TCGA cohort are shown in Fig. S2. The results demonstrated that the low-risk subtype has significantly higher specific immune cells infiltration (e.g., CD4 T cells, plasma cells, natural killer cells and dendritic cells), which was consistent with the aforementioned results.

Patients with higher TMB often have a higher number of neoantigen; therefore, these patients could obtain a stronger immune response (34,35). It was found that patients in the low-risk subgroup had significantly higher TMB compared with patients in the high-risk group (P<0.001; Fig. 4C). TMB is closely associated with genomic instability. DNA damage repair (DDR)-related genes, such as BRCA1/2, TP53 and POLE, and clinical factors were therefore considered in the multivariate logistic regression model. Following adjustment of these confounding factors, the association between the low-risk subtype and higher TMB remained notably significant (OR, 2.95; 95% CI, 1.89–4.67; P<0.001; Fig. S3). Neoantigens are generated mainly owing to non-synonymous mutations. It was found that the low-risk subtype harbored a significantly higher non-synonymous mutation load compared with the high-risk group (P<0.001; Fig. 4D).

Functional annotation of the low-risk subtype

In order to verify whether the two CRC subtypes had distinct functionalities, differential analysis and GSEA pathway analysis of genome-wide gene expression data were performed on the two CRC subtypes in TCGA CRC cohort. Analysis of differentially expressed genes in the KEGG pathway demonstrated that signaling pathways associated with immune response activation, such as ‘allograft rejection’, ‘graft versus host disease’, and ‘antigen processing and presentation’ were significantly enriched in the low-risk subtype (Fig. 5A). GO analysis also revealed the enrichment of ‘antigen processing and presentation’ signaling in the low-risk subtype. Both the ‘adaptive immune response’ and ‘immune system process’ pathways were also significantly enriched in low-risk patients (Fig. 5B). These results suggested that patients in the low-risk subgroup had a better immune response.

Figure 5. Pathway analysis of the low-risk colorectal cancer subtype, with signaling pathway annotations from (A) KEGG and (B) GO. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; NES, normalized enrichment score; FDR, false discovery rate.

Discussion

The present study identified two clinically distinct subtypes of CRC, with markedly different clinical outcomes and immune infiltration patterns, by using immune-related signatures on TCGA cohort. For responses to ICI therapy in the two CRC subtypes, the TIDE algorithm demonstrated a low score in patients of the low-risk subtype compared with that in patients of the high-risk subtype. Furthermore, higher PD-L1 expression, TIL proportion and TMB were remarkably enriched in the low-risk subtype. Finally, GSEA pathway analysis indicated immune response activation-, and ‘antigen processing and presentation’-related pathways to be significantly enriched in the low-risk subtype. Overall, the present study suggests that patients in the low-risk group have the potential to respond better to ICI treatment.

The NMF algorithm is a useful tool to perform clustering or subtyping, and extract genomic signatures; previous studies have shown its utility in distinct research directions, such as identification of image pattern, signal processing and text mining (36–38). The TIDE prediction score is a better predictor of ICI therapy compared with PD-L1 expression and TMB. The algorithm uses gene expression signatures to model two major mechanisms of tumor immune evasion: The induction of T-cell dysfunction in tumors with high infiltration of cytotoxic T lymphocytes (CTLs) and the prevention of T-cell infiltration in tumors with a low CTL level (25). In the present study, patients of the low-risk subtype harbored a significantly lower TIDE prediction score compared with patients of the high-risk subtype in TCGA and validation cohorts, suggesting that this subtype has potential for an improved response to ICI.

To date, only PD-L1 expression is confirmed by the FDA as an official criterion for ICI therapy (18). Based on the results from the KEYNOTE-001 clinical trial (18), high PD-L1 expression is now an essential condition for the use of pembrolizumab in NSCLC. The present study also demonstrated markedly higher PD-L1 expression in the low-risk CRC subtype compared with that in the high-risk subtype in TCGA and validation cohorts. TMB is emerging as a potential biomarker to predict the response to ICI therapy. Three clinical trials, including KEYNOTE-001, CHECKMATE-026 and CHECKMATE-227, demonstrated that patients with higher TMB responded more effectively to ICI (39,40). A recent study reported that high TMB was positively associated with the response to PD-1/PD-L1 inhibitors (such as pembrolizumab) in metastatic CRC, and may serve as a biomarker to predict associated immune therapy effects (41). The present study discovered that TMB in the low-risk CRC subtype was significantly higher compared with that in the high-risk subtype in TCGA cohort. Meanwhile, the non-synonymous mutation load was also markedly enriched in the low-risk CRC subgroup. Studies have reported that mutations of BRCA1/2, POLE and TP53 genes also indicate high TMB, due to the loss of DDR (42–44). Thus, multivariate analysis was conducted following the adjustment of these factors (i.e., mutations of BRCA1/2, POLE and TP53), and the association of low-risk subtype with higher TMB remained statistically significant. Several studies have revealed that TILs play crucial roles in the tumor-immune microenvironment, immune response and prognosis of CRC (44–47). In the present study, results from the CIBERSORT algorithm showed that the low-risk CRC subtype harbored a significantly higher TIL proportion compared with the high-risk group in TCGA and independent validation cohorts. Thus, TMB, TIL proportion and non-synonymous mutations load may contribute to the response to ICI in patients of the low-risk CRC subtype.

The present study has several limitations. Firstly, the gene expression data from the four validation cohorts were from different platforms, which may create bias in the analyses. Secondly, the results associated with the CRC mutation data derived from TCGA cohort were not sufficiently validated, owing to the unavailability of mutation data from the other cohorts. Thirdly, the number of CRC samples with a gene expression profile is currently insufficient to perform molecular subtyping, and thus the results could not be validated using biological experiments such as histochemical staining.

In summary, two clinically distinct CRC subtypes were identified that have significantly different survival outcomes and immune microenvironments. The low-risk CRC subtype may indicate improved response to ICI therapy, due to higher PD-L1 expression, TIL proportion and TMB. Thus, this molecular classification study and integrated multi-omic analysis of CRC may lead to a novel therapeutic approach for improving the prognosis of patients in the low-risk CRC subtype.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

XZ and GC designed the study; GC developed the methodology and acquired the related data; GC, LW, TD, YC and CC performed the data analysis and interpretation; GC and XZ drafted and revised the manuscript; XZ supervised the study. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References