The RNA sequencing (RNA-seq) gene expression data of 612 CRC samples and their clinicopathologic information were downloaded from TCGA database (https://gdc-portal.nci.nih.gov/) as our discovery cohort. RNA-seq data of 598 CRC samples (GSE39582) and the clinical information were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). An independent cohort GSE39582 was used for further validation. Building on TCGA datasets, gene mutation information was obtained via secondary databases like cBioportal (http://www.cbioportal.org/). To discover unbiased immune subtypes, we identified 1,989 IRGs in five different categories including single-cell RNA-seq data, the gene of co-stimulatory and co-inhibitory molecules, the gene of cytokine and cytokine receptors, genes involved in antigen processing and presentation, and other immune-related genes (Chen et al., 2019; Li et al., 2019).

The cBioPortal (http://www.cbioportal.org/) was used to analyze and visualize cancer genomics data sets including copy number variation and mutation (Cerami et al., 2012). The genetic alterations of potential tumor antigens against tumors were obtained from cBioPortal based on 612 CRC samples in TCGA.

The online database Gene Expression Profiling Interactive Analysis (GEPIA2) (http://gepia2.cancer-pku.cn/) which is based on samples from the TCGA and the GTEx (Genotype-Tissue Expression) databases was used for gene expression analysis to further clarify the relationship between tumor antigen gene expression and CRC prognosis (Tang et al., 2017). ANOVA was used for differential expression analysis to compare the tumor to paired normal samples. The chromosomal distribution of over- or under-expressed genes were plotted in differential genes with |Log2FC| > 1 and q-value < 0.01. Kaplan–Meier estimates of OS and RFS were performed to evaluate the prognostic value of identified tumor antigens.

To assess the correlation between the infiltration level of APCs and the expression of the identified tumor antigens TIMER database was used. (https://cistrome.shinyapps.io/timer/). TIMER includes 10,897 samples across 32 cancer types from TCGA to estimate the abundance of immune infiltrates. Prognostic gene expression data was used to calculate the abundance of six immune cell types (B cells, CD4+T cells, CD8+T cells, neutrophils, macrophages, and dendritic cells) in CRC patients using the TIMER (Li et al., 2017).

Based on the expression of 1,989 IRGs, the 201 IRGs impacting prognosis of CRC patients by univariate Cox regression analysis was performed in subsequent analysis. Consensus clustering analysis was used to identify the four immune subtypes by using the R package “Consensus Cluster Plus” with parameters defined as a maximum of nine clusters. To assess the robustness of the identified immune subtypes, we further conducted the same algorithm in the validation cohort (GEO). In-group-proportion statistic was calculated to validate the similarity of immune subtypes across the discovery and validation cohorts.

The correlation between immune subtypes and clinicopathological characteristics was evaluated by the Chi-square test. The prognostic value of the immune subtypes with OS as the endpoint was analyzed by the log-rank test and multivariable Cox proportional hazards regression model. The ANOVA algorithm was performed to assess the correlation between immune subtypes and different immune-related molecular and cellular characteristics. Finally, we employed the ssGSEA algorithm via R packages “GSVA” to comprehensively calculate the enrichment score of CRC samples on each immune-related term.

Considering the dynamic nature of the immune system, the dimensionality reduction analysis using the graph learning-based method was performed to reveal the intrinsic structure and visualize the distribution of individual patients (Trapnell et al., 2014). The discriminative dimensionality reduction with trees (DDRTree) was performed for dimension reduction (Qiu et al., 2017). Finally, the immune landscape was visualized with the function plot cell trajectory (Monocle) with the color corresponding to the immune subtype identified above (Qi et al., 2017; Wang and Mao, 2019).

All statistical analyses and data visualization were performed in R version 4.0.5, and GraphPad Prism 8.0. p < 0.05 was considered statistically significant.

To identify potential tumor antigens of CRC, 2658 over-expressed genes that could potentially encode tumor-associated antigens (TAA) were detected by screening for aberrantly expressed genes (Figure 1A). In addition, 14,751 mutated genes encoding TAA were filtered by analyzing fraction of genome altered and tumor mutational count in CRC patients (Figures 1B,C). Mutational analysis showed that the adenomatous polyposis coli (APC), tumor suppressor P53 (TP53), KRAS proto-oncogene GTPase (KARS), titin (TTN), spectrin repeat containing nuclear envelope protein 1 (SYNE1), mucin 16 (MUC16), and low-density lipoprotein receptor-related protein 1B (LRP1B) were the most frequently mutated genes in the fraction genome alteration and tumor mutational count group (Figures 1D,E). Collectively, we identified 2052 over-expressed and mutated genes.

The prognosis-associated tumor antigens were selected from the above-mentioned genes as promising candidates for designing the mRNA vaccine. Fifty-one genes were strongly correlated with the overall survival (OS) of CRC. Only six genes were closely related to the disease-free survival (RFS) (Supplementary Figure S1). The survival analysis revealed that elevated expression of thrombospondin 2 (THBS2), follistatin like 3 (FSTL3), troponin T1 (TNNT1), biglycan (BGN), collagen triple helix repeat containing 1 (CTHRC1), and NADPH oxidase 4 (NOX4) were associated with poor OS in CRC (Figures 2A–F). In addition, we found a shorter RFS time in patients with higher expression levels of THBS2, FSTL3, TNNT1, BGN, CTHRC1, and NOX4 (Supplementary Figure S2). The effectiveness of mRNA vaccines strongly relies on the availability of the target antigen(s) to professional APCs, especially dendritic cells (DCs), macrophages, and B cells. Therefore, we investigated whether the expression of THBS2, FSTL3, TNNT1, BGN, CTHRC1, and NOX4 was correlated with immune cell infiltration levels in CRC via the Tumor Immune Estimation Resource (TIMER) database. The results demonstrated that elevated expression of THBS2, FSTL3, BGN, CTHRC1, and NOX4 was exhibited in increased infiltration of macrophages, DCs and B cells in CRC, although it was insignificant in B cells. Additionally, TNNT1 expression has the same tendency of increased infiltration of macrophages, DCs, and B cells (Figures 3A–F). Taken together, we designed mRNA encoding the antigen of interest flanked by 5′ and 3′ untranslated regions (UTRs). As shown in Figure 4, six tumor antigens (THBS2, FSTL3, TNNT1, BGN, CTHRC1, and NOX4) were identified as potential therapeutic candidates for the anti-CRC mRNA vaccine that can be processed and presented by APCs to activate a robust immune response.

Immunophenotyping can be used to reflect the immune status of the tumor microenvironment, thus, allow the selection of patients for mRNA vaccine treatment. Therefore, consensus clustering was performed in 568 CRC samples based on the immune-related genes (IRGs). The results were visualized using a cumulative distribution function (CDF) plot and a delta area plot, in which k represents the number of immune subtypes. Four robust immune subtypes (IS1-IS4) were identified and the IRGs appeared to be stably clustered when K = 4 (Figures 5A–C). Significantly prognostic impact of the immune subtypes in TCGA cohort was observed. Immune subtypes IS2 and IS4 were associated with a better prognosis for OS. Conversely, patients in IS1 and IS3 displayed the poorer survival probability in the TCGA cohort (Figure 5D). Consistent with our results obtained from the TCGA cohort, each of the immune subtypes was related to prognosis in the GEO cohort as well (Figure 5E). Among all immune subtypes, IS2 and IS4 were associated with the lower stage (Figures 5F,G) compared with IS1 and IS3. Taken together, immunophenotyping could be a robust prognostic indicator and its accuracy is superior to conventional stage of the tumor in CRC.

High tumor mutational burden (TMB) can increase the tumor immunogenicity, allow the immune system to recognize and attack the cancer cells by T cell-mediated antitumor immunity. Therefore, TMB for each patient can be detected and calculated using the VarScan method in the TCGA cohort. Subtypes IS2 and IS4 had significantly higher TMB than subtypes IS1 and IS3 (Figure 6A). Furthermore, somatic mutation profiling revealed a high frequency of mutations in APC (78%), TP53 (59%), TTN (50%), KRAS (43%), SYNE1 (29%), and MUC16 (28%) across immune subtype (Figure 6B). These findings suggest that the immune subtype is a promising tool to predict TMB in CRC patients and clinical response to mRNA vaccine. Therefore, patients belong to IS2 and IS4 may respond positively to the mRNA vaccine.

Immune checkpoint molecules and human leukocyte antigen (HLA) are essential for immune function and have diverse clinical implications in cancer immunotherapy. Therefore, we next analyzed the relationship between the CRC immune subtypes and the expression level of immune checkpoint molecule genes and HLA. We found that 30 immune checkpoint-related genes were significantly differentially expressed among the immune subtypes. For example, TNFSF4, TNFRSF9, TNFRSF4, TNFRSF18, LAG3, JAK2, PTPRC, IL12B, IFNG, ICOS, HAVCR2, CTLA4, CD8A, CD86, CD80, CD40, CD28, and CD274 were significantly downregulated in IS3 tumors, while TNFSF9, TNFSF4 TNFRSF9, TNFRSF4, TNFRSF18, PTPRC, CTLA4, CD8A, CD86, CD80, CD40, CD28, and CD274 were strongly elevated in IS2 tumors (Figures 7A,B). Next, we assessed the expression level of HLA genes in the four immune subtypes and determined that the expression levels were significantly elevated in IS2 and IS4 than in IS1 and IS3 (Figures 7C,D). HLA plays important role in the activation of anti-tumor immunity through interactions with T cell receptors. Thus, the response of patients with IS2 and IS4 on mRNA vaccine treatment might be more promising.

The tumor immune microenvironment (TIME) plays a crucial role in CRC progression, thus influences the effectiveness of the mRNA vaccine. We assessed the relationship between the immune subtypes and 29 immune-related molecular features in the TCGA and GEO cohorts by the single-sample gene set enrichment analysis (ssGSEA) algorithm (Figures 8A,B). The results showed that IS2 and IS4 have significantly higher levels of B cells, CD+8 T cells, Th1 cells, and macrophages compared with IS1 and IS3. Moreover, we calculated stromal and immune scores to quantify the abundances of infiltrating stromal and immune cells for each sample in TCGA and GEO. These results also demonstrated that IS2 and IS4 have higher Immune Score than IS1 and IS3 (Figure 8C). Thus, IS2 and IS4 are immune-hot phenotypes, while IS1 and IS3 are immune-cold phenotypes. Similar results were also validated in patients in the GEO cohort (Figure 8D).

To better understand the association of CRC immune subtypes with six pan-cancer immunotypes we extracted 11 immune-related molecular features from a previous study and demonstrated that the expression of 11 molecular signatures was significantly different among the four immune subtypes (Figure 8E). Antigen-specific T cell receptor (TCR) and B cell receptor (BCR) repertoires play a significant role in the immune system for the recognition of malignant cells. We observed that IS2 and IS4 had higher scores in TCR richness and BCR richness. Higher TCR richness may enhance anticancer immunity. Therefore, immunotype can be used to reflect CRC immune status and to select suitable patients for mRNA vaccine treatment.

The immune landscape is increasingly recognized to be important in cancer development, progression, and response to therapy. Subsequently, we sought to visualize the immune landscape by the graph learning-based dimensionality reduction technique. We observed that the immune-hot subtype IS2 and immune-cold subtype IS3 were distributed at the opposite end of the horizontal axis in the immune landscape (Figure 9A). The horizontal axis in the immune landscape may represent the overall immune infiltration. As shown in Figure 9B, the horizontal axis was associated with various immune cells, of which activated CD8+ T cells, activated CD4+ T cells, and type 1 T helper cells showed the positive correlation.

The immune landscape further revealed significant intra-cluster heterogeneity within each subtype. IS3 could be further divided into four subgroups based on the distribution location in the immune landscape (Figure 9C). The further prognostic analysis demonstrated that IS3D has the best prognosis among IS3A, IS3B, and IS3C (Figure 9D). Collectively, the immune landscape based on immune subtypes can precisely recognize immune components of each CRC patient as well as predict their prognoses that are beneficial for selecting personalized therapeutics for mRNA vaccine.