1. Introduction

Colorectal cancer is the most prevalent type of malignant tumor of the gastrointestinal tract [Reference Chen, Zheng and Zong1]. In 2020, colorectal cancer is the third most common cancer, with more than 193000 diagnosed cases [Reference Xi and Xu2], leading to the second cause of cancer-related deaths, with around 830000 fatalities. Studies have found that nearly 50% of colorectal cancer patients will eventually have distant metastases, and the resistance of colorectal cancer to chemotherapy drugs has led to treatment failure and tumor metastasis invasion [Reference Filip, Vymetalkova and Petera3]. Immunotherapy, in addition to surgery and chemotherapy, has been identified as a promising treatment for specific subtypes of CRC [Reference Golshani and Zhang4]. However, there were differences in prognosis amongst patients with the same disease stage, which were linked to distinct genetic abnormalities, underlining CRC’s molecular heterogeneity.

The human immune system will go through three stages of tumor development: immune clearance, immune balance, and immune escape [Reference Gonzalez, Hagerling and Werb5]. Based on this theory, immunotherapy plays a role to restore or enhance the antitumor effect of the immune system. With the development of tumor cytology, immunotherapy has been paid more and more attention in the treatment of malignant tumors and has achieved good efficacy in some subtypes of colorectal cancer patients [Reference Ventola6]. Currently, the research related to immune checkpoint inhibitors is relatively mature and widely used in tumor immunotherapy. Immunotherapy drugs primarily target PD-1, PD-L1, and CTLA-4 [Reference Relecom, Merhi and Inchakalody7]. At present, PD-1/PD-L1 inhibitors are relatively commonly used in the immunotherapy of colorectal cancer, and their efficacy is relatively favorable in patients with deficient mismatch repair (dMMR) or microsatellite instability-high (MSI-H) types, and the efficacy in patients with other subtypes needs to be further explored, in which the level of PD-L1 expression may have an impact on the efficacy [Reference Oliveira, Bretes and Furtado8]. However, mismatch repair proficient (MMR) or microsatellite stable (MSS) types account for the majority of patients with colorectal cancer. In recent years, some researchers have begun to focus on the immune escape mechanism of pMMR/MSS colorectal cancer, trying to convert its immune “cold environment” into a “hot environment.” Recent advances in genomics and bioinformatics have facilitated identification of new immune-related genes through cancer genome sequencing [Reference Chen, Wu and Lu9]. Therefore, it is promising to develop a novel genomic and expression-based classification of CRC with an improved clinical significance.

In this study, we compared gene expression differences between wild-type genes and mutant genes using expression profiles and genomic data, as well as performed Cox regression analysis. In addition, bioinformatics analysis was employed to investigate the clinical profile and distinct characteristics of immunogenicity of different subtypes of CRC. Our findings could provide novel potential predictive indicators and therapy targets for CRC immunotherapy.

2. Materials and Methods

3. Results

3.1. The Landscape of Somatic Mutations of CRC

536 WES samples data were conducted to analyze the landscape of somatic mutations of CRC. There were 525 samples (97.95%) altered in a total of 536 samples (Figure 1(a)). In addition, 31 genes were identified as significantly mutated genes including APC, TP53, KRAS, PIK3CA, TNN, SYNE1, MUC16, FAT4, RYR2, OBSCN, ZFHX4, LRP1B, DNAH5, DNAH11, FAT3, CSMD3, FBXW7, PCLO, CSMD1, ABCA13, USH2A, RYR1, FLG, NEB, RYR3, ADGRV1, LRP2, CCDC168, DCHS2, ATM, and A1BG (Figure 1(a)). The variant classification can be divided into 9 types, among which missense mutations account for the majority (Figures 1(b) and 1(c)). The predominant somatic mutation types were C > T (Figure 1(d)). The median variants per sample were 106. TTN, APC, MUC16, SYNE1, TP53, FAT4, KARS, RYR2, OBSCN, and PIK3CA were the top 10 mutant genes (Figures 1(e) and 1(g)).

FIGURE 1: Landscape of somatic mutations of CRC. (a) Oncoplot and waterfall plot showing the somatic landscape of top 31 mutated genes in CRC. The distribution of variant classification (b), variant type (c), and SNV class (d) present. (e)-(f) Mutation load of each sample (variant classification type), and stacked bar graph shows the top ten mutated genes, including TTN, APC, MUC16, SYNE1, TP53, FAT4, KARS. RYR2, OBSCN, and PIK3CA.

3.2. Identification of Survival-Associated Somatic Mutated Genes

To investigate the link between somatic mutations and CRC prognosis, 31 significantly mutated genes were divided into two groups, respectively. We identified 4 genes SYNE1, TNN, CCDC168, and NEN mutations were significantly related to short overall survival (Figures 2(a)–2(d)). In addition, we performed the expression profiles of significantly mutated genes in one normal group and two tumor groups. We found that 20 genes were dramatically changed between the normal colon tissue and tumor group (Figure 3(a)). Among them, DNAH5, TP53, OBSCN, LRP2, NEB, PCLO, MUC16, USH2A, and CCDC168 were significantly upregulated. CSMD1, SYNE1, RYR1, RYR3, APC, ADGRV1M, DCHS2, KRAS, LRP1B, and FAT4 were significantly downregulated (Figure 3(a)) Furthermore, the relationship between the expression of somatically mutated genes and progression-free survival (PFS) in CRC was investigated. We identified high expressions of CSMD3, DCHS2, LRP2, RYR2, and ZFHX4 were significantly negatively correlated with PFS in CRC (Figures 3(b)–3(f)).

FIGURE 2: Identification of survival-associated somatic mutated genes. The CRC patients with SYNE1 (a), TNN (b), CCDC168 (c), and NEN (d) mutations had a shorter overall survival time.

FIGURE 3: The dysregulation of somatically mutated genes was correlated to the prognosis of CRC. (a) 20 genes dramatically changed between the normal colon tissue and CRC group. Among them, DNAH5, TP53, OBSCN, LRP2, NEB, PCLO, MUC16, USH2A, and CCDC168 significantly upregulated. CSMD1, SYNE1, RYR1, RYR3, APC, ADGRV1M, DCHS2, KRAS, LRP1B, and FAT4 significantly downregulated. High expression of CSMD3 (b), DCHS2 (c), LRP2 (d), RYR2 (e), and ZFHX4 (f) significantly negatively correlated with shorter PFS in CRC.

3.3. Consensus Clustering Analysis for CRC Based on the Expression of Significantly Somatic Mutated Genes

Consensus clustering was performed using the Consensus Cluster Plus R program. At consensus index, with k = 3, the cumulative distribution function has the lowest rangeability (Figure 4(a)). At k = 3, the analysis had the best delta area scores (Figure 4(b)). In total, three subtypes of CRC were identified in 620 CRC samples, including cluster1 (n = 453), cluster2 (n = 158), and cluster 3 (n = 9) based on expression level of significantly somatic mutated genes (Figures 4(c) and 4(d)). Our data show that clustered subtypes defined by the expression level of somatically mutated genes are closely related to the heterogeneity of CRC patients (Figure 5(a)).

FIGURE 4: Consensus clustering analysis for CRC based on the expression of significantly somatic mutated genes. (a) The cumulative distribution function analysis of CRC subtypes. (b) Relative change in the area under the CDF curve (CDF delta area) presented. At k = 3, the analysis had the best delta area scores. (c) Three subtypes of CRC identified in CRC samples. (d) The PCA analysis of three subtypes of CRC performed.

FIGURE 5: Clinicopathological features analysis of somatically mutated genes related subgroups of CRC. (a) Heatmap showing the expression levels of somatically mutated genes in subgroups of CRC. (b) Kaplan–Meier survival analysis of the correlation between subgroups and PFS time in CRC. (c) The distribution analysis of T stage in different subgroups of CRC. (d) The distribution analysis of N stage in different subgroups of CRC. (e) The distribution analysis of M stage in different subgroups of CRC. (f) The distribution analysis of grades in different subgroups of CRC.

3.5. Association between Immune Infiltration and Genomic Consensus Cluster in CRC

In subtyping analysis, we compared the immune infiltration of the three subgroups by the CIBERSORT algorithm. We found that B cell native, T cell CD8+, T cell CD4+ memory activated, T cell gamma delta, NK cell resting, macrophage M0/M2, myeloid dendritic cell activated, mast cell activated, and mast cell resting had differences among the three groups (Figure 6).

FIGURE 6: Association between immune infiltration and genomic consensus cluster in CRC. (a) Heatmap analysis showing immune infiltration of immune cells in three subgroups by the CIBERSORT algorithm. (b) Box-plot analysis showing immune infiltration levels in three subgroups by the CIBERSORT algorithm.

Additionally, a higher number of B cell native were produced in subtype C3 than in other subtypes; T cell CD4+ memory resting/activated and NK cell resting were significantly lower in C3 than in C1 and C2 subtypes; macrophage M0/M2, myeloid dendritic cell activated, and mast cell activated were significantly higher in C2 than in C1 and C3 subtypes; interestingly, mast cell resting was significantly lower in C2 than in C1 and C3 subtypes (Figures 6(a) and 6(b)).

4. Discussion

In the current study, we characterized the landscape of genomic alterations in CRC patients. Based on the results of WES, we identified 31 significantly mutated genes. Among them, several genes including TP53, KRAS, APC, PI3KCA, and BRAF were reported as significantly mutated genes in previous studies. In the current study, the most frequently mutated gene was TP53, which encodes tumor suppressor p53, affecting approximately 60% of CRC patients. A previous study showed that loss of p53 transcriptional activity leads to uncontrolled cell growth in various organs, including the colon [Reference Chen13]. In addition, KRAS mutations are the main intestinal cancer markers. Studies showed that 30% of human malignant tumors were related to KRAS gene mutation. The mutated KRAS is activated and unable to create normal RAS protein, causing the intracellular signal to be disrupted, resulting in uncontrolled cell proliferation and cancer [Reference Jancik, Drabek, Radzioch and Hajduch14]. Furthermore, mutations in the APC gene are the molecular pathological basis of adenomatous polyposis (FAP) and play a crucial role in the development of sporadic colorectal cancer [Reference Narayan and Roy15]. APC protein deficiency results in β-catenin accumulation in the cytoplasm; this leads to sustained transcriptional activation of TCF-mediated genes and promotes colorectal cancer progression [Reference Morin, Sparks and Korinek16]. BRAF mutations are observed in 8–12% of patients with advanced disease, and a valine amino acid substitution (V600E) in exon 15 is the most common alteration. BRAF mutations in CRC cause MAPK/ERK signaling abnormal activation, which affects cell growth and differentiation pathways [Reference Smiech, Leszczynski, Kono, Wardell and Taniguchi17]. Mutations in the PIK3CA gene cause the creation of an altered p110 subunit, which permits PI3K to hypersignal without being regulated, triggering malignant cell actions such as proliferation and migration [Reference Ligresti, Militello and Steelman18]. In addition to these commonly mutated genes, genomic mutations in SYNE1, TTN, NEB, and CCDC168 were found to be associated with a poor prognosis of colorectal cancer in this study. SYNE1 encodes a spectrin repeat-containing protein Nesprin-1 that localizes to the nuclear membrane. Mutation of SYNE1 in exonic rs9479297 leads to the upregulation of proteins SYNE1, which may contribute to cell proliferation and migration in hepatocellular and transitional cell carcinoma [Reference Chu, Kee and Lin19]. Loss of Nesprin-1 in CRC cells may alter cell destiny, contributing to carcinogenesis. In our study, somatic-mutated SYNE1 cases displayed a worse overall survival. Djulbegovic et al. demonstrated that in ocular surface squamous neoplasia, this gene mutation is the most common genetic abnormality. TTN mutations are also linked to resistance to topical interferon alpha-2b (IFN-2b) therapy, which promotes chromosomal instability, oncogenesis, and a changed response to IFN-2b treatment. Xie et al. discovered that TTN mutation was enriched in samples with high immunostimulatory signatures and that the mutation load within TTN implies high TMB status, which is consistent with our finding that somatically mutated TTN cases had a lower overall survival [Reference Xie, Tang and Sheng20]. NEB encodes nebulin, a giant protein component of the cytoskeletal matrix. In taxol-resistant ovarian cancer cells, NEB along with DCDC2, ANKRD18B, ALDH1A1, and ITGBL1 were overexpressed suggesting the involvement of these AR-related genes in taxol resistance. NEB was also identified as a potential cancer-related gene, but research into NEB in CRC is still in its early stages. In our study, we demonstrate that NEB-mutated cases are associated with the worse overall survival in CRC patients.

In addition, we performed the expression profiles of significantly mutated genes in one normal group and two tumor groups; we identified 20 DEGs; among them CSMD3, DCHS2, LRP2, RYR2, and ZFHX4 were significantly negatively correlated with PFS. Lu et al. discovered that CSMD3 is linked to tumor mutation burden and immune infiltration in ovarian cancer patients [Reference Lu, Liu and Xu21] and was identified as a potential driver gene in prostate adenocarcinoma. However, it is unclear whether CSMD3 is a driver protein in CRC and how CSMD3 contributed to tumorigenesis. DCHS2 is a big protein with numerous cadherin domains that are thought to play a role in cell adhesion. This gene may be relevant in the stomach and colorectal malignancies with high microsatellite instability, according to genome-wide association studies. It suggests that DCHS2 may relate to the immune response of MSI-H [Reference An, Je, Yoo and Lee22]. LRP2/megalin was one of the first endocytic cargos identified for the Dab2 adaptor and has been found related to poor prognosis in various cancer. In thyroid tumors, LRP2 mediated the suppressive effect of metformin on cancer proliferation by blocking JNK signaling [Reference He, Cao, Wang, Liu, Huang and Gong23]. LRP2/megalin expression influences melanoma cell growth and survival rates in frequently acquired melanoma tumors [Reference Andersen, Hammer and Hager24]. In addition, LRP2 was reported to relate to fibrosis-associated diseases and cancer through the TGF-β pathway. Furthermore, RYR2 is an important player in steroid metabolism and cancer research. Several studies uncovered that RYR2 mutation played a positive side in tumor prognosis. Wei et al. found that RYR2 3′UTR polymorphisms remained significant in the genetic susceptibility of progesterone receptor positive breast cancer [Reference Wei, Wang and Zhang25]. Liu et al. found the RYR2 mutation was linked to a greater TMB and a better clinical outcome by enhancing the antitumor immune response in breast cancer and esophageal adenocarcinoma [Reference Liu, Liu and Jiao26]. In nonsmall cell lung cancer, RYR2 mutation may prolong survival via downregulation of DKK1 and upregulation [Reference Ren, Li and Chen27]. In this study, we found that RYR2 expression was negatively related to the PFS in CRC. We speculate that the heterogeneity that exists between tumors causes RYR2 to exhibit different functions in colorectal cancer or due to the insufficient sample size of the analyzed data. ZFHX4 is a 397 kD putative transcription factor with 4 homeodomains and 22 zinc fingers that were discovered lately [Reference Chudnovsky, Kim and Zheng28]. Chudnovsky et al. reported ZFHX4 interacts with CHD4 to govern the glioblastoma tumor initiating cell state [Reference Chudnovsky, Kim and Zheng28]. ZFHX4 is also found as a prognostic factor for ovarian serous cystadenocarcinoma, esophageal cancer, and lung adenocarcinoma. Understanding the role of pivotal genetic variations in CRC cancer development and growth is a promising direction for future research.

CRC is a highly heterogeneous disease, and clinically similar tumors with similar pathology differ significantly in terms of treatment response and patient survival. Pathologic staging is critical to the efficacy of biological therapy for CRC, and the traditional classification of CRC does not fully reflect tumor heterogeneity and cannot accommodate the needs of modern cancer therapy. The concept of molecular classification is to shift the classification of tumors from morphological to molecular characteristics through comprehensive molecular analysis techniques. Afterward, with the in-depth study of the mechanism at the molecular level and the advancement of sequencing technology, the molecular typing methods of tumors have been expanded from relying on single or few markers detection to the stage of spectroscopic typing, which has led more researchers to perform molecular typing of CRC. However, there is no uniform standard for molecular typing of CRC for now. After analyzing 18 different CRC gene expression datasets, the International CRC subtyping collaborative (CRCSC) established 4 CRC molecular signature consensus subgroups (CMS) in 2015, including CMS1, CMS2, CMS3, and CMS4 [Reference Guinney, Dienstmann and Wang29]. Among them, CMS1 (MSI immune subtype) is distinguished by MSI, BRAF mutation, high CpG island methylator phenotype (CIMP), immune infiltration, and poor survival; CMS2 (typical) has a high somatic copy number alterations (SCNA) level; CMS3 (metabolic) is characterized by low SCNA and CIMP, KRAS mutation, and metabolic dysregulation; CMS4 (mesenchymal) is characterized by high SCNA, stromal infiltration, TGF activation, angiogenesis, and short survival. CMS typing is currently considered to be the most convincing method for colorectal cancer staging, and many researchers have conducted studies on CRC-targeted therapy based on this method. In this study, we identified the CRC subtype associated with the PFS. As shown in Figure 5(b), subtype C1 had the longest PFS with a median time of 8.2 years, while subtypes C2 and C3 had 4.1 and 2.7 years PFS, respectively. Patients’ clinicopathological staging was also related to the CRC subtype that we identified.

Following the determination of the three subtypes based on the expression of CRC-related mutated genes, the immune infiltration of the three groups was evaluated by using the CIBERSORT method. We demonstrated that B cell native, T cell CD8+, T cell CD4+ memory activated, T cell gamma delta, NK cell resting, macrophage M0, macrophage M2, myeloid dendritic cell activated, mast cell activated, and mast cell resting had differences among the three groups. According to the findings, the three subgroups classified by 22 somatically significantly mutated genes had a high capacity to differentiate patients with different immune statuses, which is helpful for the prediction of immunotherapy response of CRC patients.

This study still has some limitations. First, despite the data from CRC patients used in this study being extracted from TCGA, the sample size was still small; second, despite the obvious clinical significance of the newly proposed subtypes, little is known about their underlying mechanism. As a result, we must investigate the molecular mechanisms of the three subtypes. Third, there is no validation of the functions of hub genes.

Finally, the current investigation confirmed the predictive significance of somatically altered genes and offered a novel genomic categorization with clinical relevance. This discovery presented a foundation for CRC research and molecular classification of CRC types to guide precision therapy by interpreting genomic data.